v4.3

delicatessen/estimating_equations/__init__.py

@@ -13,7 +13,7 @@
                           )
 
 from .regression import (ee_regression, ee_glm, ee_mlogit, ee_beta_regression, ee_tobit,
-                         ee_robust_regression,
+                         ee_robust_regression, ee_expectile_regression,
                          ee_ridge_regression, ee_lasso_regression, ee_dlasso_regression,
                          ee_elasticnet_regression, ee_bridge_regression,
                          ee_additive_regression,

delicatessen/estimating_equations/regression.py

@@ -777,6 +777,135 @@ def ee_tobit(theta, X, y, lower=None, upper=None, weights=None, offset=None):
     return np.vstack([ef_treg, ef_sigma])
 
 
+def ee_expectile_regression(theta, X, y, model, tau=0.5, weights=None, offset=None):
+    r"""Estimating equations for expectile regression. Expectile regression is a generalization of linear (and
+    non-linear) mean regression models for other parts of the distributions. These other parts are referred to as
+    'expectiles', which are analogous to quantiles (but not the same). Expectile regression is implemented through
+    an asymmetric weight, :math:`\tau` on the negative and positive residuals. The general estimating equation is
+
+    .. math::
+
+        \sum_{i=1}^n \omega_i \times \left\{ Y_i - g(X_i^T \theta) \right\} X_i = 0
+
+    where :math:`g` indicates a transformation function. Here, :math:`g` is the identity function for linear,
+    inverse-logit function, :math:`\text{expit}(u) = 1 / (1 + \exp(u))`, for logistic, and :math:`\exp(u)` for Poisson.
+    The parameters are :math:`\theta`, which is a 1-by-`b` array corresponding to the coefficients in the regression
+    model and `b` is the distinct covariates included as part of ``X``. For example, if ``X`` is a 3-by-`n` matrix, then
+    :math:`\theta` will be a 1-by-3 array.
+
+    Here, :math:`\omega_i` are the asymmetric weights, where a weight of :math:`\tau` is assigned when the residual of
+    an observation is above 0 and :math:`1 - \tau` otherwise. Note that :math:`0 < \tau < 1` for expectile regression.
+    While related to quantiles (and thus quantile regression), there are distinctions. Expectile regression is less
+    straightforward to interpret, but it is computationally convenient to compute due to it retaining the L2 loss
+    function. See the references for further details.
+
+    Note
+    ----
+    At :math:`\tau = 0.5`, expectile regression is equivalent to regression for the mean.
+
+    Parameters
+    ----------
+    theta : ndarray, list, vector
+        Theta in this case consists of `b` values. Therefore, initial values should consist of the same number as the
+        number of columns present. This can easily be implemented by ``[0, ] * X.shape[1]``.
+    X : ndarray, list, vector
+        2-dimensional vector of `n` observed values for `b` variables.
+    y : ndarray, list, vector
+        1-dimensional vector of `n` observed values.
+    model : str
+        Type of regression model to estimate. Options are ``'linear'`` (linear regression), ``'logistic'`` (logistic
+        regression), and ``'poisson'`` (Poisson regression).
+    tau : float, optional
+        Asymmetric weight for expectile regression. This parameter controls which of the expectiles is being estimated.
+        Default is ``0.5``, which defaults to linear regression. Note that :math:`0 < \tau < 1`.
+    weights : ndarray, list, vector, None, optional
+        1-dimensional vector of `n` weights. Default is ``None``, which assigns a weight of 1 to all observations.
+    offset : ndarray, list, vector, None, optional
+        A 1-dimensional offset to be included in the model. Default is ``None``, which applies no offset term.
+
+    Examples
+    --------
+    Construction of an estimating equation(s) with ``ee_expectile_regression`` should be done similar to the following
+
+    >>> import numpy as np
+    >>> import pandas as pd
+    >>> from scipy.stats import logistic
+    >>> from delicatessen import MEstimator
+    >>> from delicatessen.estimating_equations import ee_expectile_regression
+
+    Some generic data to estimate the regression models
+
+    >>> n = 500
+    >>> data = pd.DataFrame()
+    >>> data['X'] = np.random.normal(size=n)
+    >>> data['Z'] = np.random.normal(size=n)
+    >>> data['Y'] = 0.5 + 2*data['X'] - 1*data['Z'] + np.random.normal(loc=0, size=n)
+    >>> data['C'] = 1
+
+    Note that ``C`` here is set to all 1's. This will be the intercept in the regression.
+
+    >>> X_mat = data[['C', 'X', 'Z']]
+    >>> y = data['Y']
+
+    To start, we will demonstrate expectile regression for the continuous outcome ``Y``. Defining psi, or the stacked
+    estimating equations
+
+    >>> def psi(theta):
+    >>>     return ee_expectile_regression(theta=theta, X=X_mat, y=y, model='linear', tau=0.5)
+
+    Calling the M-estimator (note that ``init`` requires 3 values, since ``X.shape[1]`` is 3).
+
+    >>> estr = MEstimator(stacked_equations=psi, init=[0., 0., 0.,])
+    >>> estr.estimate()
+
+    Inspecting the parameter estimates, variance, and confidence intervals
+
+    >>> estr.theta
+    >>> estr.variance
+    >>> estr.confidence_intervals()
+
+    This previous specification was the same as linear regression. Now we will estimate the 0.75 expectile
+
+    >>> def psi(theta):
+    >>>         return ee_expectile_regression(theta=theta, X=X_mat, y=y, model='linear', tau=0.75)
+
+    >>> estr = MEstimator(stacked_equations=psi, init=[0., 0., 0.,])
+    >>> estr.estimate()
+
+    Weighted models can be estimated by specifying the optional ``weights`` argument.
+
+    References
+    ----------
+    Newey WK, & Powell JL. (1987). Asymmetric least squares estimation and testing. *Econometrica*, 55(4), 819-847.
+
+    Waltrup LS, Sobotka F, Kneib T, & Kauermann G. (2015). Expectile and quantile regression—David and Goliath?.
+    *Statistical Modelling*, 15(5), 433-456.
+    """
+    # Input checking for expectile requested
+    if not 0 < tau < 1:
+        raise ValueError("Expectile regression is only defined for values of tau within the interval (0,1). The "
+                         "specified tau, " + str(tau) + ", lies outside this interval.")
+
+    # Preparation of input shapes and object types
+    X, y, beta, offset = _prep_inputs_(X=X, y=y, theta=theta, penalty=None, offset=offset)
+
+    # Applying transformation function for corresponding model
+    transform = _model_transform_(model=model)    # Looking up corresponding transformation
+    pred_y = transform(np.dot(X, beta) + offset)  # Generating predicted values via speedy matrix calculation
+    residual = y - pred_y                         # Residual from regression model
+
+    # Asymmetric weight for corresponding expectile
+    asym_w = np.where(residual >= 0, tau, 1 - tau)
+
+    # Apply optional user-specified weights
+    if weights is not None:
+        w = generate_weights(weights=weights, n_obs=X.shape[0])
+        asym_w = asym_w * w[:, None]
+
+    # Output b-by-n matrix
+    return (asym_w * residual * X).T              # Return weighted expectile regression score function
+
+
 #################################################################
 # Robust Regression Estimating Equations
 

delicatessen/estimation.py

@@ -48,6 +48,7 @@ def __init__(self, stacked_equations, init, subset=None, finite_correction=None)
         self._dx_ = ''            # Derivative approximation error
         self._pinv_ = ''          # Pseudo-inverse
         self._over_identified_ = None
+        self._theta_buffer_ = np.asarray(self.init, dtype=float)
 
     def confidence_intervals(self, alpha=0.05):
         r"""Calculate two-sided Wald-type :math:`(1 - \alpha) \times` 100% confidence intervals using the point
@@ -412,23 +413,15 @@ def _eval_ee_(stacked_equations, subset):
         numpy.array
         """
         # IF stacked_equation returns 1 value, only return that 1 value
-        if len(stacked_equations.shape) == 1:          # Checking length
-            vals = np.sum(stacked_equations)           # ... avoid SciPy error by returning value rather than tuple
+        if len(stacked_equations.shape) == 1:                # Checking length of estimating function parameters
+            vals = np.sum(stacked_equations)                 # ... avoid SciPy error by return value rather than tuple
         # ELSE need to return a tuple for the root-finding procedure
-        else:                                          # ... also considering how subset argument is handled
-            # NOTE: switching to np.sum(..., axis=1) didn't speed things up versus a for-loop
-            vals = ()                                  # ... create empty tuple
-            rows = stacked_equations.shape[0]          # ... determine how many rows / parameters are present
-            if subset is None:                         # ... if no subset, then simple loop where
-                for i in range(rows):                  # ... go through each individual theta in the stack
-                    row = stacked_equations[i, :]      # ... extract corresponding row
-                    vals += (np.sum(row), )            # ... then add the theta sum to the tuple of thetas
-            else:                                      # ... if subset, then conditional loop (to speed up)
-                for i in range(rows):                  # ... go through each individual theta in the stack
-                    if i in subset:                    # ... if parameter is in subset then
-                        row = stacked_equations[i, :]  # ... extract corresponding row
-                        vals += (np.sum(row), )        # ... then add the theta sum to the tuple of thetas
-            vals = np.asarray(vals)                    # ... converting to a NumPy array for ease
+        else:                                                # ... also considering how subset argument is handled
+            if subset is None:                               # ... if no subset,
+                efuncs = stacked_equations                   # ... then keep all estimating functions
+            else:                                            # ... if subset,
+                efuncs = stacked_equations[list(subset), :]  # ... then keep only the subset of estimating functions
+            vals = np.sum(efuncs, axis=1)                    # ... Transform EFs to EEs
 
         # Return the calculated values of theta
         return vals
@@ -780,14 +773,16 @@ def _mestimation_answer_(self, theta):
         if self._subset_ is None:                      # If NOT subset then,
             full_theta = theta                         # ... then use the full input theta
         else:                                          # If subset then,
-            full_theta = np.asarray(self.init)         # ... copy the initial values to ndarray
-            np.put(full_theta,                         # ... update in place the previous array
-                   ind=self._subset_,                  # ... go to the subset indices
-                   v=theta)                            # ... then input current iteration values
+            # Then in _mestimation_answer_:
+            self._theta_buffer_[:] = self.init         # ... carry over intial values to buffer array
+            np.put(self._theta_buffer_,                # ... update in place the previous array
+                   ind=self._subset_,                  # ... at the subset indices
+                   v=theta)                            # ... with current iteration solver values
+            full_theta = self._theta_buffer_           # ... update name convention to match no-subset
 
         stacked_equations = np.asarray(self.stacked_equations(full_theta))  # Returning stacked equation
-        return self._eval_ee_(stacked_equations=stacked_equations,   # Passing to evaluating function
-                              subset=self._subset_)                  # ... with specified subset
+        return self._eval_ee_(stacked_equations=stacked_equations,          # Passing to evaluating function
+                              subset=self._subset_)                         # ... with specified subset
 
     @staticmethod
     def _solve_coefficients_(stacked_equations, init, method, maxiter, tolerance):
@@ -1219,10 +1214,11 @@ def _gmmestimation_answer_(self, theta):
         if self._subset_ is None:                      # If NOT subset then,
             full_theta = theta                         # ... then use the full input theta
         else:                                          # If subset then,
-            full_theta = np.asarray(self.theta)        # ... copy the initial values to ndarray
-            np.put(full_theta,                         # ... update in place the previous array
-                   ind=self._subset_,                  # ... go to the subset indices
-                   v=theta)                            # ... then input current iteration values
+            self._theta_buffer_[:] = self.init         # ... carry over intial values to buffer array
+            np.put(self._theta_buffer_,                # ... update in place the previous array
+                   ind=self._subset_,                  # ... at the subset indices
+                   v=theta)                            # ... with current iteration solver values
+            full_theta = self._theta_buffer_           # ... update name convention to match no-subset
 
         # Evaluating estimating functions
         stacked_equations = np.asarray(self.stacked_equations(full_theta))  # Returning stacked equation

delicatessen/version.py

@@ -1 +1 @@
-__version__ = "4.2"
+__version__ = "4.3"