Repository: pzivich/Delicatessen (language: python)

Issue:
Tests for writing own equations, built-in equations, one-liner lambda functions



Discussion / hints:

Fix patch:
diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml
index b1854d1..1644c79 100644
--- a/.github/workflows/python-package.yml
+++ b/.github/workflows/python-package.yml
@@ -27,7 +27,7 @@ jobs:
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
-        python -m pip install pytest numpy scipy statsmodels
+        python -m pip install pytest numpy scipy statsmodels scikit-learn
     - name: Install package
       run: |
         pip install -e .
diff --git a/.gitignore b/.gitignore
index 7e7c1d2..c90541c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,4 +3,5 @@ __pycache__/
 builddir/
 build/
 dist/
+.idea/
 
diff --git a/README.md b/README.md
index cf0aa6e..a4a0c8e 100644
--- a/README.md
+++ b/README.md
@@ -2,9 +2,10 @@
 
 # Delicatessen
 
+![tests](https://github.com/pzivich/Delicatessen/actions/workflows/python-package.yml/badge.svg)
 [![version](https://badge.fury.io/py/delicatessen.svg)](https://badge.fury.io/py/delicatessen)
 [![docs](https://readthedocs.org/projects/deli/badge/?version=latest)](https://deli.readthedocs.io/en/latest/?badge=latest)
-![tests](https://github.com/pzivich/Delicatessen/actions/workflows/python-package.yml/badge.svg)
+[![Downloads](https://pepy.tech/badge/delicatessen/month)](https://pepy.tech/project/delicatessen)
 
 The one-stop sandwich (variance) shop in Python. `delicatessen` is a Python 3.6+ library for the generalized calculus 
 of M-estimation. 
@@ -13,44 +14,46 @@ of M-estimation.
 but the releases may not be currently stable until the true release. I would not recommend on integrating
 `delicatessen` into your full workflow at this time.
 
+## M-Estimation and Estimating Equations
+
 Here, we provide a brief overview of M-estimation theory. For a more detailed and formal introduction to M-estimation,
 I highly recommend chapter 7 of Boos & Stefanski (2013). M-estimation is a generalization of robust inference (here
 robust refers to allowing for misspecification of secondary assumptions does not invalidate inference) for
 likelihood-based methods to a general context. *M-estimators* are solutions to estimating equations. A large number of 
 consistent and asymptotically normal statistics can be put into the M-Estimation framework. Some examples include: 
-mean, regression, delta method, and many others.
+mean, regression, delta method, and among others.
 
 To apply the M-Estimator, we solve the stacked estimating equations using observed data. This is similar to other 
-approaches, but the key advantage of M-Estimators is the straightforward estimation of the variance following from this 
-framework, under suitable regularity conditions. Specifically, the sandwich variance.
+approaches, but the key advantage of M-Estimators is the straightforward estimation of the variance via the sandwich 
+variance.
 
 While M-Estimation is a powerful tool, the derivatives and matrix algebra can quickly become unwieldy. This is where 
-`delicatessen` comes in. `delicatessen` takes stacked estimating equations and data and works through all the necessary 
-calculations. Therefore, M-Estimation can be more widely adopted without needing to solve every derivative for your 
-particular problem. We can let the computer do all that hard work of finding the roots, numerically approximating 
-the derivatives, and doing the matrix calculations for us.
+`delicatessen` comes in. `delicatessen` takes an array of stacked estimating equations and data and works through the 
+root-finding, numerically approximating the partial derivatives, and matrix calculations. Therefore, M-Estimation can 
+be more widely adopted without needing to solve every derivative for your particular problem. We can let the computer 
+do all that hard math for us.
 
 In addition to implementing a general M-estimator, `delicatessen` also comes with a variety of built-in estimating 
 equations. See the [delicatessen website](https://deli.readthedocs.io/en/latest/) for the full set of available
 estimating equations and how to use them.
 
-# Installation
+## Installation
 
-## Installing:
+### Installing:
 
 You can install via `python -m pip install delicateseen`
 
-## Dependencies:
+### Dependencies:
 
 There are only two dependencies: `numpy`, `scipy`
 
 To replicate the tests located in `tests/`, you will additionally need to install: `panda`, `statsmodels`, and `pytest`
 
-While `delicatessen` will likely work with older versions of NumPy and SciPy (this has not been formally tested),
-it is recommended to use `numpy >= 1.18.0` and `scipy >= 1.4.0`. There has been no reported testing on versions
-earlier than these.
+While `delicatessen` is expected to work with older versions of NumPy and SciPy, this has not been formally tested.
+Therefore, it is recommended to use `numpy >= 1.18.0` and `scipy >= 1.4.0` as there is no currently reported testing 
+on previous versions.
 
-# Getting started
+## Getting started
 
 To demonstrate `delicatessen`, below is a simple demonstration of calculating the mean for the following data
 
@@ -59,7 +62,8 @@ import numpy as np
 y = np.array([1, 2, 3, 1, 4, 1, 3, -2, 0, 2])
 ```
 
-Loading the M-estimator functionality from deli, building the estimating equation, and 
+Loading the M-estimator functionality from deli, building the estimating equation, and printing the results to the
+console
 
 ```python
 from delicatessen import MEstimator
@@ -74,5 +78,5 @@ print(mestimate.theta)     # Estimate of the mean
 print(mestimate.variance)  # Variance estimator for the mean
 ```
 
-For recommendations and general advice on using `delicatessen`, see the full documentation and worked examples 
-available at [delicatessen website](https://deli.readthedocs.io/en/latest/).
+For full details on using `delicatessen`, see the full documentation and worked examples available 
+at [delicatessen website](https://deli.readthedocs.io/en/latest/).
diff --git a/delicatessen/data.py b/delicatessen/data.py
new file mode 100644
index 0000000..085f5fe
--- /dev/null
+++ b/delicatessen/data.py
@@ -0,0 +1,104 @@
+import numpy as np
+
+
+def load_shaq_free_throws():
+    """Load example data from Boos and Stefanski (2013) on Shaquille O'Neal free throws in the 2000 NBA playoffs
+    (Table 7.1 on pg 324).
+
+    Notes
+    -----
+    From left to right, the columns in the array correspond to:
+        * game - game number
+        * ft_success - free throws made during game
+        * ft_attempt - free throws attempted during game
+
+    Returns
+    -------
+    ndarray
+    """
+    d = np.array([[ 1,  4,  5],
+                  [ 2,  5, 11],
+                  [ 3,  5, 14],
+                  [ 4,  5, 12],
+                  [ 5,  2,  7],
+                  [ 6,  7, 10],
+                  [ 7,  6, 14],
+                  [ 8,  9, 15],
+                  [ 9,  4, 12],
+                  [10,  1,  4],
+                  [11, 13, 27],
+                  [12,  5, 17],
+                  [13,  6, 12],
+                  [14,  9,  9],
+                  [15,  7, 12],
+                  [16,  3, 10],
+                  [17,  8, 12],
+                  [18,  1,  6],
+                  [19, 18, 39],
+                  [20,  3, 13],
+                  [21, 10, 17],
+                  [22,  1,  6],
+                  [23,  3, 12], ])
+    return d
+
+
+def load_inderjit():
+    """Load example data from Inderjit et al. (2002) on the dose-response of herbicide on perennial ryegrass growth
+
+    Notes
+    -----
+    From left to right, the columns in the array correspond to:
+        * response - ryegrass root length
+        * dose - herbicide dose
+
+    Returns
+    -------
+    ndarray
+    """
+    d = np.array([[7.5800000,  0.00],
+                  [8.0000000,  0.00],
+                  [8.3285714,  0.00],
+                  [7.2500000,  0.00],
+                  [7.3750000,  0.00],
+                  [7.9625000,  0.00],
+                  [8.3555556,  0.94],
+                  [6.9142857,  0.94],
+                  [7.7500000,  0.94],
+                  [6.8714286,  1.88],
+                  [6.4500000,  1.88],
+                  [5.9222222,  1.88],
+                  [1.9250000,  3.75],
+                  [2.8857143,  3.75],
+                  [4.2333333,  3.75],
+                  [1.1875000,  7.50],
+                  [0.8571429,  7.50],
+                  [1.0571429,  7.50],
+                  [0.6875000, 15.00],
+                  [0.5250000, 15.00],
+                  [0.8250000, 15.00],
+                  [0.2500000, 30.00],
+                  [0.2200000, 30.00],
+                  [0.4400000, 30.00], ])
+    return d
+
+
+def load_robust_regress(outlier=True):
+    """Load illustrative example for robust linear regression.
+
+    Parameters
+    ----------
+
+    Returns
+    -------
+
+    """
+    height = [168.519, 166.944, 164.327, 164.058, 166.212, 167.358,
+              165.244, 169.352, 159.386, 166.953, 163.876,
+              164.245, 165.061, 162.876, 164.185]
+    weight = [67.634, 67.418, 63.394, 66.18, 65.573, 67.66, 64.592,
+              68.4, 62.043, 67.093, 65.202, 64.328, 64.754,
+              62.179, 64.716]
+    if outlier:
+        weight[8] += 3
+
+    return np.array([height, weight]).T
diff --git a/delicatessen/estimating_equations.py b/delicatessen/estimating_equations.py
index 943b8b3..37c30c8 100644
--- a/delicatessen/estimating_equations.py
+++ b/delicatessen/estimating_equations.py
@@ -1,5 +1,5 @@
 import numpy as np
-from delicatessen.utilities import logit, inverse_logit
+from delicatessen.utilities import logit, inverse_logit, identity
 
 #################################################################
 # Basic Estimating Equations
@@ -15,15 +15,14 @@ def ee_mean(theta, y):
     Note
     ----
     All provided estimating equations are meant to be wrapped inside a user-specified function. Throughtout, these
-    user-defined functions are defined as `psi`.
-
+    user-defined functions are defined as ``psi``.
 
     Parameters
     ----------
-    theta : vector
+    theta : ndarray, list, vector
         Theta in the case of the mean consists of a single value. Therefore, an initial value like the form of
         [0, ] is recommended.
-    y : vector
+    y : ndarray, list, vector
         1-dimensional vector of n observed values. No missing data should be included (missing data may cause unexpected
         behavior when attempting to calculate the mean).
 
@@ -34,7 +33,7 @@ def ee_mean(theta, y):
 
     Examples
     --------
-    Construction of a estimating equation(s) with `ee_mean` should be done similar to the following
+    Construction of a estimating equation(s) with ``ee_mean`` should be done similar to the following
 
     >>> from delicatessen import MEstimator
     >>> from delicatessen.estimating_equations import ee_mean
@@ -50,52 +49,53 @@ def ee_mean(theta, y):
 
     Calling the M-estimation procedure
 
-    >>> mestimation = MEstimator(stacked_equations=psi, init=[0, ])
-    >>> mestimation.estimate()
+    >>> estr = MEstimator(stacked_equations=psi, init=[0, ])
+    >>> estr.estimate()
 
     Inspecting the parameter estimates, the variance, and the asymptotic variance
 
-    >>> mestimation.theta
-    >>> mestimation.variance
-    >>> mestimation.asymptotic_variance
+    >>> estr.theta
+    >>> estr.variance
+    >>> estr.asymptotic_variance
 
     References
     ----------
     Boos DD, & Stefanski LA. (2013). M-estimation (estimating equations). In Essential Statistical Inference
     (pp. 297-337). Springer, New York, NY.
     """
-    # Output 1-by-n array
-    return np.asarray(y) - theta     # Estimating equation for the mean
+    # Convert input y values to NumPy array
+    y_array = np.asarray(y)
+
+    # Output 1-by-n array estimating equation for the mean
+    return y_array - theta
 
 
 def ee_mean_robust(theta, y, k):
-    r""" Default stacked estimating equation for robust mean (robust location) estimator. The estimating equation for
-    the robust mean is
+    r""" Default stacked estimating equation for robust mean (location) estimator. The estimating equation for the
+    robust mean is
 
     .. math::
 
         \sum_i^n \psi(Y_i, \theta_1) = \sum_i^n Y^*_i - \theta_1 = 0
 
-    where Y* is bounded between k and -k.
+    where :math:`Y^*` is bounded between :math:`k` and :math:`-k`.
 
     Note
     ----
-    Since psi is non-differentiable at k or -k, it must be assumed that the mean is sufficiently far from k. Otherwise,
-    difficulties might arise in the variance calculation.
-
+    Since psi is non-differentiable at :math:`k` or :math:`-k`, it must be assumed that the mean is sufficiently far
+    from :math:`k`. Otherwise, difficulties might arise in the variance calculation.
 
     Note
     ----
     All provided estimating equations are meant to be wrapped inside a user-specified function. Throughtout, these
-    user-defined functions are defined as `psi`.
-
+    user-defined functions are defined as ``psi``.
 
     Parameters
     ----------
-    theta : vector
+    theta : ndarray, list, vector
         Theta in the case of the robust mean consists of a single value. Therefore, an initial value like the form of
-        [0, ] is recommended.
-    y : vector
+        ``[0, ]`` is recommended.
+    y : ndarray, vector, list
         1-dimensional vector of n observed values. No missing data should be included (missing data may cause unexpected
         behavior when attempting to calculate the robust mean).
     k : int, float
@@ -108,7 +108,7 @@ def ee_mean_robust(theta, y, k):
 
     Examples
     --------
-    Construction of a estimating equation(s) with `ee_mean_robust` should be done similar to the following
+    Construction of a estimating equation(s) with ``ee_mean_robust`` should be done similar to the following
 
     >>> from delicatessen import MEstimator
     >>> from delicatessen.estimating_equations import ee_mean_robust
@@ -124,14 +124,14 @@ def ee_mean_robust(theta, y, k):
 
     Calling the M-estimation procedure
 
-    >>> mestimation = MEstimator(stacked_equations=psi, init=[0, ])
-    >>> mestimation.estimate()
+    >>> estr = MEstimator(stacked_equations=psi, init=[0, ])
+    >>> estr.estimate()
 
     Inspecting the parameter estimates, the variance, and the asymptotic variance
 
-    >>> mestimation.theta
-    >>> mestimation.variance
-    >>> mestimation.asymptotic_variance
+    >>> estr.theta
+    >>> estr.variance
+    >>> estr.asymptotic_variance
 
     References
     ----------
@@ -141,310 +141,1200 @@ def ee_mean_robust(theta, y, k):
     Huber PJ. (1992). Robust estimation of a location parameter. In Breakthroughs in statistics (pp. 492-518).
     Springer, New York, NY.
     """
-    var = np.asarray(y)                   # Convert y to NumPy array
-    var = np.where(var > k, k, var)       # Apply the upper bound
-    var = np.where(var < -k, -k, var)     # Apply the lower bound
+    # Convert input y values to NumPy array
+    y_array = np.asarray(y)
+
+    # Bounding via np.clip
+    y_bound = np.clip(y_array, a_min=-k, a_max=k)
+
+    # Output 1-by-n array estimating equation for robust mean
+    return y_bound - theta
+
+
+def ee_mean_variance(theta, y):
+    r"""Default stacked estimating equation for mean and variance. The estimating equations for the mean and
+     variance are
+
+    .. math::
+
+        \sum_i^n \psi_1(Y_i, \theta_1) = \sum_i^n Y_i - \theta_1 = 0
+
+        \sum_i^n \psi_2(Y_i, \theta_1) = \sum_i^n (Y_i - \theta_1)^2 - \theta_2 = 0
+
+    Unlike ``ee_mean``, theta consists of 2 elements. The output covariance matrix will also provide estimates for each
+    of the theta values.
+
+    Note
+    ----
+    All provided estimating equations are meant to be wrapped inside a user-specified function. Throughtout, these
+    user-defined functions are defined as ``psi``.
+
+
+    Parameters
+    ----------
+    theta : ndarray, list, vector
+        Theta in this case consists of two values. Therefore, initial values like the form of [0, 0] is recommended.
+    y : ndarray, list, vector
+        1-dimensional vector of n observed values. No missing data should be included (missing data may cause unexpected
+        behavior when attempting to calculate the mean).
+
+    Returns
+    -------
+    array :
+        Returns a 2-by-n NumPy array evaluated for the input theta and y
+
+    Examples
+    --------
+    Construction of a estimating equation(s) with ``ee_mean_variance`` should be done similar to the following
+
+    >>> from delicatessen import MEstimator
+    >>> from delicatessen.estimating_equations import ee_mean_variance
+
+    Some generic data to estimate the mean for
+
+    >>> y_dat = [1, 2, 4, 1, 2, 3, 1, 5, 2]
+
+    Defining psi, or the stacked estimating equations
+
+    >>> def psi(theta):
+    >>>     return ee_mean_variance(theta=theta, y=y_dat)
+
+    Calling the M-estimation procedure (note that `init` has 2 values now).
+
+    >>> estr = MEstimator(stacked_equations=psi, init=[0, 0, ])
+    >>> estr.estimate()
+
+    Inspecting the parameter estimates, the variance, and the asymptotic variance
+
+    >>> estr.theta
+    >>> estr.variance
+    >>> estr.asymptotic_variance
+
+    For this estimating equation, ``mestimation.theta[1]`` and ``mestimation.asymptotic_variance[0][0]`` are expected
+    to always be equal.
+
+    References
+    ----------
+    Boos DD, & Stefanski LA. (2013). M-estimation (estimating equations). In Essential Statistical Inference
+    (pp. 297-337). Springer, New York, NY.
+    """
+    # Convert input y values to NumPy array
+    y_array = np.asarray(y)
+
+    # Output 2-by-n matrix of estimating equations
+    return (y_array - theta[0],                  # Estimating equation for mean
+            (y_array - theta[0])**2 - theta[1])  # Estimating equation for variance
+
+
+def ee_percentile(theta, y, q):
+    r"""Default stacked estimating equation for percentiles (or quantiles).
+
+    Note
+    ----
+    Due to this estimating equation being non-smooth, estimated percentile values may differ from the closed-form
+    definition of the percentile. In general, closed form solutions for percentiles will be preferred, but this
+    estimating equation is offered for completeness.
+
+    .. math::
+
+        \sum_i^n \psi_q(Y_i, \theta_q) = \sum_i^n q - I(Y_i \le \theta_q) = 0
+
+    Notice that this estimating equation is non-smooth. Therefore, optimization and numerically approximating
+    derivatives for this estimating equation are more difficult.
+
+    Note
+    ----
+    The following optional parameters ``MEstimator.estimate()`` may benefit from these changes ``solver='hybr'``,
+    ``dx=1``, ``order=15``, and increasing the ``tolerance``.
+
+    Parameters
+    ----------
+    theta : ndarray, list, vector
+        Theta in this case consists of two values. Therefore, initial values like the form of [0, 0] is recommended.
+    y : ndarray, list, vector
+        1-dimensional vector of n observed values. No missing data should be included (missing data may cause unexpected
+        behavior when attempting to calculate the mean).
+    q : float
+        Percentile to calculate. Must be (0, 1)
+
+    Returns
+    -------
+    array :
+        Returns a 1-by-n NumPy array evaluated for the input theta and y
+
+    Examples
+    --------
+    Construction of a estimating equation(s) with ``ee_percentile`` should be done similar to the following
+
+    >>> from delicatessen import MEstimator
+    >>> from delicatessen.estimating_equations import ee_percentile
+
+    Some generic data to estimate the mean for
+
+    >>> np.random.seed(89041)
+    >>> y_dat = np.random.normal(size=100)
+
+    Defining psi, or the stacked estimating equations
+
+    >>> def psi(theta):
+    >>>     return ee_percentile(theta=theta, y=y_dat, q=0.5)
+
+    Calling the M-estimation procedure (note that ``init`` has 2 values now).
+
+    >>> estr = MEstimator(stacked_equations=psi, init=[0, ])
+    >>> estr.estimate(solver='hybr', tolerance=1e-3, dx=1, order=15)
+
+    Notice that we use a different solver, tolerance values, and parameters for numerically approximating the derivative
+    here. These changes generally work better for percentile optimizations since the estimating equation is non-smooth.
+    Furthermore, optimization is hard when only a few observations (<100) are available. In general, closed form
+    solutions for percentiles will be preferred.
+
+    >>> estr.theta
+
+    Then displays the estimated percentile / median. In this example, there is a difference between the closed form
+    solution (``-0.07978``) and M-Estimation (``-0.06022``). Again, this results from the non-smooth estimating
+    equation.
+
+    References
+    ----------
+    Boos DD, & Stefanski LA. (2013). M-estimation (estimating equations). In Essential Statistical Inference
+    (pp. 297-337). Springer, New York, NY.
+    """
+    if q >= 1 or q <= 0:
+        raise ValueError("`q` must be (0, 1)")
+
+    # Convert input y values to NumPy array
+    y_array = np.asarray(y)
+
+    # Output 1-by-n array of the estimating equations
+    return q - 1*(y_array <= theta)
+
+
+def ee_positive_mean_deviation(theta, y):
+    """Default stacked estimating equations for the positive mean deviation. The estimating equations are
+
+    .. math::
+
+        \sum_i^n \psi_1(Y_i, \theta) = \sum_i^n 2(Y_i - \theta_2)I(Y_i > \theta_2) - \theta_1 = 0
+
+        \sum_i^n \psi_2(Y_i, \theta) = \sum_i^n 0.5 - I(Y_i \le - \theta_2) = 0
+
+    where the first estimating equation is for the positive mean difference, and the second estimating equation is for
+    the median.
+
+    Note
+    ----
+    All provided estimating equations are meant to be wrapped inside a user-specified function. Throughtout, these
+    user-defined functions are defined as ``psi``.
+
+    Parameters
+    ----------
+    theta : ndarray, list, vector
+        Theta in this case consists of two values. Therefore, initial values like the form of ``[0, 0]`` are
+        recommended.
+    y : ndarray, list, vector
+        1-dimensional vector of n observed values. No missing data should be included (missing data may cause unexpected
+        behavior when attempting to calculate the positive mean deviation).
+
+    Returns
+    -------
+    array :
+        Returns a 2-by-n NumPy array evaluated for the input theta and y
+
+    Examples
+    --------
+    Construction of a estimating equation(s) with ``ee_positive_mean_deviation`` should be done similar to the following
+
+    >>> from delicatessen import MEstimator
+    >>> from delicatessen.estimating_equations import ee_positive_mean_deviation
+
+    Some generic data to estimate the mean for
+
+    >>> y_dat = [1, 2, 4, 1, 2, 3, 1, 5, 2]
+
+    Defining psi, or the stacked estimating equations
+
+    >>> def psi(theta):
+    >>>     return ee_positive_mean_deviation(theta=theta, y=y_dat)
+
+    Calling the M-estimation procedure (note that ``init`` has 2 values now).
+
+    >>> estr = MEstimator(stacked_equations=psi, init=[0, 0, ])
+    >>> estr.estimate()
+
+    Inspecting the parameter estimates, the variance, and the asymptotic variance
+
+    >>> estr.theta
+    >>> estr.variance
+    >>> estr.asymptotic_variance
+
+    References
+    ----------
+    Boos DD, & Stefanski LA. (2013). M-estimation (estimating equations). In Essential Statistical Inference
+    (pp. 297-337). Springer, New York, NY.
+    """
+    # Convert input y values to NumPy array
+    y_array = np.asarray(y)
+
+    # Calculating median with built-in estimating equation
+    median = ee_percentile(theta=theta[1], y=y_array, q=0.5)
+
+    # Output 2-by-n matrix of estimating equations
+    return ((2*(y_array - theta[1])*(y_array > theta[1])) - theta[0],   # Estimating equation for positive mean dev
+            median, )                                                   # Estimating equation for median
+
+
+#################################################################
+# Regression Estimating Equations
+
+
+def ee_linear_regression(theta, X, y, weights=None):
+    r"""Default stacked estimating equation for linear regression without the homoscedastic assumption. The estimating
+    equation is
+
+    .. math::
+
+        \sum_i^n \psi(Y_i, X_i, \theta) = \sum_i^n (Y_i - X_i^T \theta) X_i = 0
+
+    Here, theta is a 1-by-b array, where b is the distinct covariates included as part of X. For example, if X is a
+    3-by-n matrix, then theta will be a 1-by-3 array. The code is general to allow for an arbitrary number of X's (as
+    long as there is enough support in the data).
+
+    Note
+    ----
+    All provided estimating equations are meant to be wrapped inside a user-specified function. Throughtout, these
+    user-defined functions are defined as ``psi``.
+
+    Here, theta corresponds to the coefficients in a linear regression model
+
+    Note
+    ----
+    For complex regression problems, the optimizer behind the scenes is not particularly robust (unlike functions
+    specializing in solely OLS). Therefore, optimization of OLS via a separate functionality can be done then those
+    estimated parameters are fed forward as the initial values (which should result in a more stable optimization).
+
+    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 accomplished generally by ``[0, ] * X.shape[1]``.
+    X : ndarray, list, vector
+        2-dimensional vector of n observed values for b variables. No missing data should be included (missing data
+        may cause unexpected behavior).
+    y : ndarray, list, vector
+        1-dimensional vector of n observed values. No missing data should be included (missing data may cause unexpected
+        behavior).
+    weights : ndarray, list, vector, None, optional
+        1-dimensional vector of n weights. No missing weights should be included. Default is None, which assigns a
+        weight of 1 to all observations.
+
+    Returns
+    -------
+    array :
+        Returns a b-by-n NumPy array evaluated for the input theta and y
+
+    Examples
+    --------
+    Construction of a estimating equation(s) with ``ee_linear_regression`` should be done similar to the following
+
+    >>> import numpy as np
+    >>> import pandas as pd
+    >>> from delicatessen import MEstimator
+    >>> from delicatessen.estimating_equations import ee_linear_regression
+
+    Some generic data to estimate a linear regresion model
+
+    >>> 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.
+
+    Defining psi, or the stacked estimating equations
+
+    >>> def psi(theta):
+    >>>         return ee_linear_regression(theta=theta, X=data[['C', 'X', 'Z']], y=data['Y'])
+
+    Calling the M-estimation procedure (note that ``init`` has 3 values now, since ``X.shape[1] = 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()
+
+    References
+    ----------
+    Boos DD, & Stefanski LA. (2013). M-estimation (estimating equations). In Essential Statistical Inference
+    (pp. 297-337). Springer, New York, NY.
+    """
+    X = np.asarray(X)                       # Convert to NumPy array
+    y = np.asarray(y)[:, None]              # Convert to NumPy array and ensure correct shape for matrix algebra
+    beta = np.asarray(theta)[:, None]       # Convert to NumPy array and ensure correct shape for matrix algebra
+
+    # Allowing for a weighted linear model
+    if weights is None:                     # If weights is unspecified
+        w = np.ones(X.shape[0])                 # ... assign weight of 1 to all observations
+    else:                                   # Otherwise
+        w = np.asarray(weights)                 # ... set weights as input vector
+
+    # Output b-by-n matrix
+    return w*((y -                   # Speedy matrix algebra for regression
+               np.dot(X, beta))      # ... linear regression requires no transformations
+              * X).T                 # ... multiply by coefficient and transpose for correct orientation
+
+
+def ee_robust_linear_regression(theta, X, y, k, weights=None):
+    """Default stacked estimating equation for robust linear regression. Specifically, robust linear regression is
+    robust to outlying observations of the outcome variable (``y``). The estimating equation is
+
+    .. math::
+
+        \sum_i^n \psi(Y_i, X_i, \theta) = \sum_i^n \psi_k(Y_i - X_i^T \theta) X_i = 0
+
+    where k indicates the upper and lower bounds. Here, theta is a 1-by-b array, where b is the distinct covariates
+    included as part of X. For example, if X is a 3-by-n matrix, then theta will be a 1-by-3 array. The code is general
+    to allow for an arbitrary number of X's (as long as there is enough support in the data).
+
+    Note
+    ----
+    All provided estimating equations are meant to be wrapped inside a user-specified function. Throughtout, these
+    user-defined functions are defined as ``psi``.
+
+    Here, theta corresponds to the coefficients in a robust linear regression model
+
+    Note
+    ----
+    For complex regression problems, the optimizer behind the scenes is not particularly robust (unlike functions
+    specializing in solely OLS). Therefore, optimization of OLS via a separate functionality can be done then those
+    estimated parameters are fed forward as the initial values (which should result in a more stable optimization).
+
+    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 accomplished generally by ``[0, ] * X.shape[1]``.
+    X : ndarray, list, vector
+        2-dimensional vector of n observed values for b variables. No missing data should be included (missing data
+        may cause unexpected behavior).
+    y : ndarray, list, vector
+        1-dimensional vector of n observed values. No missing data should be included (missing data may cause unexpected
+        behavior).
+    k : int, float
+        Value to set the symmetric maximum upper and lower bounds on the difference between the observations and
+        predicted values
+    weights : ndarray, list, vector, None, optional
+        1-dimensional vector of n weights. No missing weights should be included. Default is None, which assigns a
+        weight of 1 to all observations.
+
+    Returns
+    -------
+    array :
+        Returns a b-by-n NumPy array evaluated for the input theta and y
+
+    Examples
+    --------
+    Construction of a estimating equation(s) with ``ee_robust_linear_regression`` should be done similar to the
+    following
+
+    >>> import numpy as np
+    >>> import pandas as pd
+    >>> from delicatessen import MEstimator
+    >>> from delicatessen.estimating_equations import ee_robust_linear_regression
+
+    Some generic data to estimate a robust linear regresion model
+
+    >>> n = 100
+    >>> 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, scale=3, size=n)
+    >>> data['C'] = 1
+
+    Note that ``C`` here is set to all 1's. This will be the intercept in the regression.
+
+    Defining psi, or the stacked estimating equations
+
+    >>> def psi(theta):
+    >>>         return ee_robust_linear_regression(theta=theta, X=data[['C', 'X', 'Z']], y=data['Y'], k=3)
+
+
+    Calling the M-estimation procedure (note that ``init`` has 3 values now, since ``X.shape[1] = 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()
+
+    References
+    ----------
+    Boos DD, & Stefanski LA. (2013). M-estimation (estimating equations). In Essential Statistical Inference
+    (pp. 297-337). Springer, New York, NY.
+    """
+    X = np.asarray(X)                       # Convert to NumPy array
+    y = np.asarray(y)[:, None]              # Convert to NumPy array and ensure correct shape for matrix algebra
+    beta = np.asarray(theta)[:, None]       # Convert to NumPy array and ensure correct shape for matrix algebra
+
+    # Allowing for a weighted linear model
+    if weights is None:                     # If weights is unspecified
+        w = np.ones(X.shape[0])                 # ... assign weight of 1 to all observations
+    else:                                   # Otherwise
+        w = np.asarray(weights)                 # ... set weights as input vector
+
+    # Generating predictions and applying Huber function for robust
+    preds = np.clip(y - np.dot(X, beta), -k, k)
+
+    # Output b-by-n matrix
+    return w*(preds                  # ... linear regression requires no transformations
+              * X).T                 # ... multiply by coefficient and transpose for correct orientation
+
+
+def ee_ridge_linear_regression(theta, y, X, penalty=1.0, weights=None):
+    r"""Default stacked estimating equation for ridge linear regression. Ridge regression applies an L2-regularization
+    through a squared magnitude penalty. The estimating equation is
+
+    .. math::
+
+        \sum_i^n \psi(Y_i, X_i, \theta) = \sum_i^n (Y_i - X_i^T \theta) X_i - \lambda \theta = 0
+
+    Here, theta is a 1-by-b array, where b is the distinct covariates included as part of X. For example, if X is a
+    3-by-n matrix, then theta will be a 1-by-3 array. The code is general to allow for an arbitrary number of X's (as
+    long as there is enough support in the data).
+
+    The 'strength' of the penalty term is indicated by :math:`\lambda`, which is the ``penalty`` argument scaled (or
+    divided by) the number of observations.
+
+    Note
+    ----
+    All provided estimating equations are meant to be wrapped inside a user-specified function. Throughtout, these
+    user-defined functions are defined as ``psi``.
+
+    Here, theta corresponds to the coefficients in a linear regression model
+
+    Note
+    ----
+    For complex regression problems, the optimizer behind the scenes is not particularly robust (unlike functions
+    specializing in solely OLS). Therefore, optimization of OLS via a separate functionality can be done then those
+    estimated parameters are fed forward as the initial values (which should result in a more stable optimization).
+
+    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 accomplished generally by ``[0, ] * X.shape[1]``.
+    X : ndarray, list, vector
+        2-dimensional vector of n observed values for b variables. No missing data should be included (missing data
+        may cause unexpected behavior).
+    y : ndarray, list, vector
+        1-dimensional vector of n observed values. No missing data should be included (missing data may cause unexpected
+        behavior).
+    penalty : int, float, optional
+        Penalty parameter for the ridge regression, which is further scaled by the number of observations.
+        Default is 1.0
+    weights : ndarray, list, vector, None, optional
+        1-dimensional vector of n weights. No missing weights should be included. Default is None, which assigns a
+        weight of 1 to all observations.
+
+    Returns
+    -------
+    array :
+        Returns a b-by-n NumPy array evaluated for the input theta and y
+
+    Examples
+    --------
+    Construction of a estimating equation(s) with ``ee_ridge_linear_regression`` should be done similar to the
+    following
+
+    >>> import numpy as np
+    >>> import pandas as pd
+    >>> from delicatessen import MEstimator
+    >>> from delicatessen.estimating_equations import ee_ridge_linear_regression
+
+    Some generic data to estimate a linear regresion model
+
+    >>> 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.
+
+    Defining psi, or the stacked estimating equations
+
+    >>> def psi(theta):
+    >>>         return ee_ridge_linear_regression(theta=theta, X=data[['C', 'X', 'Z']], y=data['Y'], penalty=5.5)
+
+    Calling the M-estimation procedure (note that ``init`` has 3 values now, since ``X.shape[1] = 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()
+
+    References
+    ----------
+    ...
+    """
+    X = np.asarray(X)                       # Convert to NumPy array
+    y = np.asarray(y)[:, None]              # Convert to NumPy array and ensure correct shape for matrix algebra
+    beta = np.asarray(theta)[:, None]       # Convert to NumPy array and ensure correct shape for matrix algebra
+
+    # Allowing for a weighted linear model
+    if weights is None:                     # If weights is unspecified
+        w = np.ones(X.shape[0])                 # ... assign weight of 1 to all observations
+    else:                                   # Otherwise
+        w = np.asarray(weights)                 # ... set weights as input vector
+
+    # Creating penalty term for ridge regression
+    penalty_scaled = penalty / X.shape[0]         # Scaling penalties by 1/N
+    penalty_terms = (penalty_scaled *             # Penalty term applied to regression coefficients
+                     np.asarray(theta))[:, None]
+
+    # Output b-by-n matrix
+    return (w*((y -                  # Speedy matrix algebra for regression
+               np.dot(X, beta))      # ... linear regression requires no transformations
+               * X).T                # ... multiply by coefficient and transpose for correct orientation
+            - penalty_terms)         # Subtract off penalty term(s) from each observation
+
+
+def ee_logistic_regression(theta, X, y, weights=None):
+    r"""Default stacked estimating equation for logistic regression. The estimating equation is
+
+    .. math::
+
+        \sum_i^n \psi(Y_i, X_i, \theta) = \sum_i^n (Y_i - expit(X_i^T \theta)) X_i = 0
+
+    where expit, or the inverse logit is
+
+    .. math::
+
+        expit(u) = 1 / (1 + exp(u))
+
+    Here, theta is a 1-by-b array, where b is the distinct covariates included as part of X. For example, if X is a
+    3-by-n matrix, then theta will be a 1-by-3 array. The code is general to allow for an arbitrary number of X's (as
+    long as there is enough support in the data).
+
+    Note
+    ----
+    All provided estimating equations are meant to be wrapped inside a user-specified function. Throughtout, these
+    user-defined functions are defined as ``psi``.
+
+    Here, theta corresponds to the coefficients in a logistic regression model, and therefore are the log-odds.
+
+    Note
+    ----
+    For complex regression problems, the optimizer behind the scenes is not particularly robust (unlike functions
+    specializing in solely logistic regression). Therefore, optimization of logistic regression via a separate
+    functionality can be done then those estimated parameters are fed forward as the initial values (which should
+    result in a more stable optimization).
+
+
+    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 accomplished generally by ``[0, ] * X.shape[1]``.
+    X : ndarray, list, vector
+        2-dimensional vector of n observed values for b variables. No missing data should be included (missing data
+        may cause unexpected behavior).
+    y : ndarray, list, vector
+        1-dimensional vector of n observed values. The Y values should all be 0 or 1. No missing data should be
+        included (missing data may cause unexpected behavior).
+    weights : ndarray, list, vector, None, optional
+        1-dimensional vector of n weights. No missing weights should be included. Default is None, which assigns a
+        weight of 1 to all observations.
+
+    Returns
+    -------
+    array :
+        Returns a b-by-n NumPy array evaluated for the input theta and y
+
+    Examples
+    --------
+    Construction of a estimating equation(s) with ``ee_logistic_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_logistic_regression
+
+    Some generic data to estimate a logistic regresion model
+
+    >>> n = 500
+    >>> data = pd.DataFrame()
+    >>> data['X'] = np.random.normal(size=n)
+    >>> data['Z'] = np.random.normal(size=n)
+    >>> data['Y'] = np.random.binomial(n=1, p=logistic.cdf(0.5 + 2*data['X'] - 1*data['Z']), size=n)
+    >>> data['C'] = 1
+
+    Note that ``C`` here is set to all 1's. This will be the intercept in the regression.
+
+    Defining psi, or the stacked estimating equations
+
+    >>> def psi(theta):
+    >>>         return ee_logistic_regression(theta=theta, X=data[['C', 'X', 'Z']], y=data['Y'])
+
+
+    Calling the M-estimation procedure (note that `init` has 3 values now, since ``X.shape[1] = 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()
+
+    References
+    ----------
+    Boos DD, & Stefanski LA. (2013). M-estimation (estimating equations). In Essential Statistical Inference
+    (pp. 297-337). Springer, New York, NY.
+    """
+    X = np.asarray(X)                      # Convert to NumPy array
+    y = np.asarray(y)[:, None]             # Convert to NumPy array and ensure correct shape for matrix algebra
+    beta = np.asarray(theta)[:, None]      # Convert to NumPy array and ensure correct shape for matrix algebra
+
+    # Allowing for a weighted linear model
+    if weights is None:                     # If weights is unspecified
+        w = np.ones(X.shape[0])                 # ... assign weight of 1 to all observations
+    else:                                   # Otherwise
+        w = np.asarray(weights)                 # ... set weights as input vector
+
+    # Output b-by-n matrix
+    return w*((y -                                # Speedy matrix algebra for regression
+               inverse_logit(np.dot(X, beta)))    # ... inverse-logit transformation of predictions
+              * X).T                              # ... multiply by coefficient and transpose for correct orientation
+
+
+#################################################################
+# Dose-Response Estimating Equations
+
+
+def ee_4p_logistic(theta, X, y):
+    r"""Default stacked estimating equation estimating equations for the four parameter logistic model (4PL). 4PL is
+    often used for dose-response and bioassay analyses. The estimating equations are
+
+    .. math::
+
+        \sum_i^n \psi(Y_i, X_i, \theta) = \sum_i^n (Y_i - expit(X_i^T \theta)) X_i = 0
+
+        \sum_i^n \psi(Y_i, X_i, \theta) = \sum_i^n (Y_i - expit(X_i^T \theta)) X_i = 0
+
+        \sum_i^n \psi(Y_i, X_i, \theta) = \sum_i^n (Y_i - expit(X_i^T \theta)) X_i = 0
+
+        \sum_i^n \psi(Y_i, X_i, \theta) = \sum_i^n (Y_i - expit(X_i^T \theta)) X_i = 0
+
+    Here, theta is a 1-by-4 array, where 4 are the 4 parameters of the 4PL. The first theta corresponds to lower limit,
+    the second corresponds to the effective dose (ED50), the third corresponds to the steepness of the curve, and the
+    fourth corresponds to the upper limit.
+
+    Note
+    ----
+    All provided estimating equations are meant to be wrapped inside a user-specified function. Throughtout, these
+    user-defined functions are defined as ``psi``.
+
+    Parameters
+    ----------
+    theta : ndarray, list, vector
+        Theta in this case consists of 4 values. In general, starting values ``>0`` are better choices for the 4PL model
+    X : ndarray, list, vector
+        1-dimensional vector of n dose values. No missing data should be included (missing data may cause unexpected
+        behavior).
+    y : ndarray, list, vector
+        1-dimensional vector of n response values. No missing data should be included (missing data may cause
+        unexpected behavior).
+
+    Returns
+    -------
+    array :
+        Returns a 4-by-n NumPy array evaluated for the input theta, y, X
+
+    Examples
+    --------
+    Construction of a estimating equation(s) with ``ee_4p_logistic`` should be done similar to the following
+
+    >>> from delicatessen import MEstimator
+    >>> from delicatessen.data import load_inderjit
+    >>> from delicatessen.estimating_equations import ee_4p_logistic
+
+    For demonstration, we use dose-response data from Inderjit et al. (2002), which can be loaded from ``delicatessen``
+    directly.
+
+    >>> d = load_inderjit()   # Loading array of data
+    >>> dose_data = d[:, 1]   # Dose data
+    >>> resp_data = d[:, 0]   # Response data
+
+    Defining psi, or the stacked estimating equations
+
+    >>> def psi(theta):
+    >>>     return ee_4p_logistic(theta=theta, X=dose_data, y=resp_data)
+
+    The 4PL model and others are harder to optimize compared to other estimating equations. Namely, the optimizer is
+    not aware of implicit bounds on the parameters. To reduce non-convergence issues, we can give the root-finder good
+    starting values.
+
+    For the 4PL, the upper limit should *always* be greater than the lower limit. Second, the ED50 should be between
+    the lower and upper limits. Third, the sign for the steepness depends on whether the response declines (positive)
+    or the response increases (negative). Finally, some solvers may be better suited to the problem, so try a few
+    different options.
+
+    Here, we use some general starting values that should perform well in many cases. For the lower-bound, give the
+    minimum response value as the initial. For ED50, give the mid-point between the maximum response and the minimum
+    response. The initial value for steepness is more difficult. Ideally, we would give a starting value of zero, but
+    that will fail in this example. Giving a small positive starting value works in this example. For the upper-bound,
+    give the maximum response value as the initial. Finally, we use the ``lm`` solver.
+
+    Note
+    ----
+    To summarize the recommendations, be sure to examine your data (e.g., scatterplot). This will help to determine the
+    initial starting values for the root-finding procedure. Otherwise, you may come across a convergence error.
+
+    >>> estr = MEstimator(psi, init=[np.min(resp_data),
+    >>>                              (np.max(resp_data)+np.min(resp_data)) / 2,
+    >>>                              (np.max(resp_data)+np.min(resp_data)) / 2,
+    >>>                              np.max(resp_data)])
+    >>> estr.estimate(solver='lm')
+
+    Inspecting the parameter estimates, variance, and confidence intervals
+
+    >>> estr.theta
+    >>> estr.variance
+    >>> estr.confidence_intervals()
+
+    Inspecting the parameter estimates
+
+    >>> estr.theta[0]    # lower limit
+    >>> estr.theta[1]    # ED(50)
+    >>> estr.theta[2]    # steepness
+    >>> estr.theta[3]    # upper limit
+
+    References
+    ----------
+    Ritz C, Baty F, Streibig JC, & Gerhard D. (2015). Dose-response analysis using R. *PloS One*, 10(12), e0146021.
+
+    An H, Justin TL, Aubrey GB, Marron JS, & Dittmer DP. (2019). dr4pl: A Stable Convergence Algorithm for the 4
+    Parameter Logistic Model. *R J.*, 11(2), 171.
+
+    Inderjit, Streibig JC, & Olofsdotter M. (2002). Joint action of phenolic acid mixtures and its significance in
+    allelopathy research. *Physiologia Plantarum*, 114(3), 422-428.
+    """
+    # Creating rho to cut down on typing
+    rho = (X / theta[1]) ** theta[2]
+
+    # Generalized 4PL model function for y-hat
+    fx = theta[0] + (theta[3] - theta[0]) / (1 + rho)
+
+    # Using a special implementatin of natural log here
+    nested_log = np.log(X / theta[1],             # ... to avoid dose=0 issues only take log
+                        where=0 < X)              # ... where dose>0 (otherwise puts zero in place)
 
-    # Output 1-by-n array
-    return var - theta                    # Estimating equation for robust mean
+    # Calculate the derivatives for the gradient
+    deriv = np.array((1 - 1/(1+rho),                                           # Gradient for lower limit
+                     (theta[3]-theta[0])*theta[2]/theta[1]*rho/(1+rho)**2,     # Gradient for steepness
+                     (theta[3] - theta[0]) * nested_log * rho / (1 + rho)**2,  # Gradient for ED50
+                     1 / (1 + rho)), )                                         # Gradient for upper limit
 
+    # Compute gradient and return for each i
+    return -2*(y-fx)*deriv
 
-def ee_mean_variance(theta, y):
-    r"""Default stacked estimating equation for mean and variance. The estimating equations for the mean and
-     variance are
+
+def ee_3p_logistic(theta, X, y, lower):
+    r"""Default stacked estimating equation estimating equations for the three parameter logistic model (3PL). 3PL is
+    often used for dose-response and bioassay analyses. The estimating equations are
 
     .. math::
 
-        \sum_i^n \psi_1(Y_i, \theta_1) = \sum_i^n Y_i - \theta_1 = 0
+        \sum_i^n \psi(Y_i, X_i, \theta) = \sum_i^n (Y_i - expit(X_i^T \theta)) X_i = 0
 
-        \sum_i^n \psi_2(Y_i, \theta_1) = \sum_i^n (Y_i - \theta_1)^2 - \theta_2 = 0
+        \sum_i^n \psi(Y_i, X_i, \theta) = \sum_i^n (Y_i - expit(X_i^T \theta)) X_i = 0
 
-    Unlike `ee_mean`, theta consists of 2 elements. The output covariance matrix will also provide estimates for each
-    of the theta values.
+        \sum_i^n \psi(Y_i, X_i, \theta) = \sum_i^n (Y_i - expit(X_i^T \theta)) X_i = 0
+
+    Here, theta is a 1-by-3 array, where 3 are the 3 parameters of the 3PL. The first theta corresponds to the
+    effective dose (ED50), the second corresponds to the steepness of the curve, and the third corresponds to the upper
+    limit. The lower limit is pre-specified by the user (and is no longer estimated)
 
     Note
     ----
     All provided estimating equations are meant to be wrapped inside a user-specified function. Throughtout, these
-    user-defined functions are defined as `psi`.
-
+    user-defined functions are defined as ``psi``.
 
     Parameters
     ----------
-    theta : vector
-        Theta in this case consists of two values. Therefore, initial values like the form of [0, 0] is recommended.
-    y : vector
-        1-dimensional vector of n observed values. No missing data should be included (missing data may cause unexpected
-        behavior when attempting to calculate the mean).
+    theta : ndarray, list, vector
+        Theta in this case consists of 3 values. In general, starting values ``>0`` are better choices for the 3PL model
+    X : ndarray, list, vector
+        1-dimensional vector of n dose values. No missing data should be included (missing data may cause unexpected
+        behavior).
+    y : ndarray, list, vector
+        1-dimensional vector of n response values. No missing data should be included (missing data may cause
+        unexpected behavior).
+    lower : int, float
+        Set value for the lower limit.
 
     Returns
     -------
     array :
-        Returns a 2-by-n NumPy array evaluated for the input theta and y
+        Returns a 3-by-n NumPy array evaluated for the input theta, y, X
 
     Examples
     --------
-    Construction of a estimating equation(s) with `ee_mean_variance` should be done similar to the following
+    Construction of a estimating equation(s) with ``ee_3p_logistic`` should be done similar to the following
 
     >>> from delicatessen import MEstimator
-    >>> from delicatessen.estimating_equations import ee_mean_variance
+    >>> from delicatessen.data import load_inderjit
+    >>> from delicatessen.estimating_equations import ee_3p_logistic
 
-    Some generic data to estimate the mean for
+    For demonstration, we use dose-response data from Inderjit et al. (2002), which can be loaded from ``delicatessen``
+    directly.
 
-    >>> y_dat = [1, 2, 4, 1, 2, 3, 1, 5, 2]
+    >>> d = load_inderjit()   # Loading array of data
+    >>> dose_data = d[:, 1]   # Dose data
+    >>> resp_data = d[:, 0]   # Response data
 
-    Defining psi, or the stacked estimating equations
+    Since there is a natural lower-bound of 0 for root growth, we set ``lower=0``. Defining psi, or the stacked
+    estimating equations
 
     >>> def psi(theta):
-    >>>     return ee_mean_variance(theta=theta, y=y_dat)
+    >>>     return ee_3p_logistic(theta=theta, X=dose_data, y=resp_data,
+    >>>                           lower=0)
 
-    Calling the M-estimation procedure (note that `init` has 2 values now).
+    The 3PL model and others are harder to optimize compared to other estimating equations. Namely, the optimizer is
+    not aware of implicit bounds on the parameters. To reduce non-convergence issues, we can give the root-finder good
+    starting values.
 
-    >>> mestimation = MEstimator(stacked_equations=psi, init=[0, 0, ])
-    >>> mestimation.estimate()
+    For the 3PL, the upper limit should *always* be greater than the set lower limit. Second, the ED50 should be between
+    the lower and upper limits. Third, the sign for the steepness depends on whether the response declines (positive)
+    or the response increases (negative). Finally, some solvers may be better suited to the problem, so try a few
+    different options.
 
-    Inspecting the parameter estimates, the variance, and the asymptotic variance
+    Here, we use some general starting values that should perform well in many cases. For ED50, give the mid-point
+    between the maximum response and the minimum response. The initial value for steepness is more difficult. Ideally,
+    we would give a starting value of zero, but that will fail in this example. Giving a small positive starting value
+    works in this example. For the upper-bound, give the maximum response value as the initial. Finally, we use the
+    ``lm`` solver.
 
-    >>> mestimation.theta
-    >>> mestimation.variance
-    >>> mestimation.asymptotic_variance
+    Note
+    ----
+    To summarize the recommendations, be sure to examine your data (e.g., scatterplot). This will help to determine the
+    initial starting values for the root-finding procedure. Otherwise, you may come across a convergence error.
+
+    >>> estr = MEstimator(psi, init=[(np.max(resp_data)+np.min(resp_data)) / 2,
+    >>>                              (np.max(resp_data)+np.min(resp_data)) / 2,
+    >>>                              np.max(resp_data)])
+    >>> estr.estimate(solver='lm')
+
+    Inspecting the parameter estimates, variance, and confidence intervals
 
-    In this example, `mestimation.theta[1]` and `mestimation.asymptotic_variance[0][0]` are expected to be equal.
+    >>> estr.theta
+    >>> estr.variance
+    >>> estr.confidence_intervals()
+
+    Inspecting the parameter estimates
+
+    >>> estr.theta[0]    # ED(50)
+    >>> estr.theta[1]    # steepness
+    >>> estr.theta[2]    # upper limit
 
     References
     ----------
-    Boos DD, & Stefanski LA. (2013). M-estimation (estimating equations). In Essential Statistical Inference
-    (pp. 297-337). Springer, New York, NY.
+    Ritz C, Baty F, Streibig JC, & Gerhard D. (2015). Dose-response analysis using R. *PloS One*, 10(12), e0146021.
+
+    An H, Justin TL, Aubrey GB, Marron JS, & Dittmer DP. (2019). dr4pl: A Stable Convergence Algorithm for the 4
+    Parameter Logistic Model. *R J.*, 11(2), 171.
+
+    Inderjit, Streibig JC, & Olofsdotter M. (2002). Joint action of phenolic acid mixtures and its significance in
+    allelopathy research. *Physiologia Plantarum*, 114(3), 422-428.
     """
-    # Output 2-by-n matrix
-    return (y - theta[0],                  # Estimating equation for mean
-            (y - theta[0])**2 - theta[1])  # Estimating equation for variance
+    # Creating rho to cut down on typing
+    rho = (X / theta[0])**theta[1]
 
+    # Generalized 3PL model function for y-hat
+    fx = lower + (theta[2] - lower) / (1 + rho)
 
-#################################################################
-# Regression Estimating Equations
+    # Using a special implementation of natural log here
+    nested_log = np.log(X / theta[0],             # ... to avoid dose=0 issues only take log
+                        where=0 < X)              # ... where dose>0 (otherwise puts zero in place)
 
+    # Calculate the derivatives for the gradient
+    deriv = np.array(((theta[2]-lower)*theta[1]/theta[0]*rho/(1+rho)**2,     # Gradient for steepness
+                      (theta[2]-lower) * nested_log * rho / (1+rho)**2,      # Gradient for ED50
+                      1 / (1 + rho)), )                                      # Gradient for upper limit
 
-def ee_linear_regression(theta, X, y, weights=None):
-    r"""Default stacked estimating equation for linear regression. The estimating equation is
+    # Compute gradient and return for each i
+    return -2*(y - fx)*deriv
 
-    .. math::
 
-        \sum_i^n \psi(Y_i, X_i, \theta) = \sum_i^n (Y_i - X_i^T \theta) X_i = 0
+def ee_2p_logistic(theta, X, y, lower, upper):
+    r"""Default stacked estimating equation estimating equations for the two parameter logistic model (2PL). 2PL is
+    often used for dose-response and bioassay analyses. The estimating equations are
 
-    Here, theta is a 1-by-b array, where b is the distinct covariates included as part of X. For example, if X is a
-    3-by-n matrix, then theta will be a 1-by-3 array. The code is general to allow for an arbitrary number of X's (as
-    long as there is enough support in the data).
+    .. math::
 
-    Note
-    ----
-    All provided estimating equations are meant to be wrapped inside a user-specified function. Throughtout, these
-    user-defined functions are defined as `psi`.
+        \sum_i^n \psi(Y_i, X_i, \theta) = \sum_i^n (Y_i - expit(X_i^T \theta)) X_i = 0
+
+        \sum_i^n \psi(Y_i, X_i, \theta) = \sum_i^n (Y_i - expit(X_i^T \theta)) X_i = 0
 
+        \sum_i^n \psi(Y_i, X_i, \theta) = \sum_i^n (Y_i - expit(X_i^T \theta)) X_i = 0
 
-    Here, theta corresponds to the coefficients in a linear regression model
+    Here, theta is a 1-by-2 array, where 2 are the 2 parameters of the 2PL. The first theta corresponds to the
+    effective dose (ED50), and the second corresponds to the steepness of the curve. Both the lower limit and upper
+    limit are pre-specified by the user (and no longer estimated).
 
     Note
     ----
-    For complex regression problems, the optimizer behind the scenes is not particularly robust (unlike functions
-    specializing in solely OLS). Therefore, optimization of OLS via a separate functionality can be done then those
-    estimated parameters are fed forward as the initial values (which should result in a more stable optimization).
+    All provided estimating equations are meant to be wrapped inside a user-specified function. Throughtout, these
+    user-defined functions are defined as ``psi``.
 
     Parameters
     ----------
-    theta : 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 accomplished generally by `[0, ] * X.shape[1]`.
-    X : vector
-        2-dimensional vector of n observed values for b variables. No missing data should be included (missing data
-        may cause unexpected behavior).
-    y : vector
-        1-dimensional vector of n observed values. No missing data should be included (missing data may cause unexpected
+    theta : ndarray, list, vector
+        Theta in this case consists of 2 values. In general, starting values >0 are better choices for the 3PL model
+    X : ndarray, list, vector
+        1-dimensional vector of n dose values. No missing data should be included (missing data may cause unexpected
         behavior).
-    weights : vector, None, optional
-        1-dimensional vector of n weights. No missing weights should be included. Default is None, which assigns a
-        weight of 1 to all observations.
+    y : ndarray, list, vector
+        1-dimensional vector of n response values. No missing data should be included (missing data may cause
+        unexpected behavior).
+    lower : int, float
+        Set value for the lower limit.
+    upper : int, float
+        Set value for the upper limit.
 
     Returns
     -------
     array :
-        Returns a b-by-n NumPy array evaluated for the input theta and y
+        Returns a 2-by-n NumPy array evaluated for the input theta, y, X
 
     Examples
     --------
-    Construction of a estimating equation(s) with `ee_linear_regression` should be done similar to the following
+    Construction of a estimating equation(s) with ``ee_2p_logistic`` should be done similar to the following
 
-    >>> import numpy as np
-    >>> import pandas as pd
     >>> from delicatessen import MEstimator
-    >>> from delicatessen.estimating_equations import ee_linear_regression
-
-    Some generic data to estimate a linear regresion model
+    >>> from delicatessen.data import load_inderjit
+    >>> from delicatessen.estimating_equations import ee_2p_logistic
 
-    >>> 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
+    For demonstration, we use dose-response data from Inderjit et al. (2002), which can be loaded from ``delicatessen``
+    directly.
 
-    Note that `C` here is set to all 1's. This will be the intercept in the regression.
+    >>> d = load_inderjit()   # Loading array of data
+    >>> dose_data = d[:, 1]   # Dose data
+    >>> resp_data = d[:, 0]   # Response data
 
-    Defining psi, or the stacked estimating equations
+    Since there is a natural lower-bound of 0 for root growth, we set ``lower=0``. While a natural upper bound does not
+    exist for this example, we set ``upper=8`` for illustrative purposes. Defining psi, or the stacked estimating
+    equations
 
     >>> def psi(theta):
-    >>>         return ee_linear_regression(theta=theta, X=data[['C', 'X', 'Z']], y=data['Y'])
+    >>>     return ee_2p_logistic(theta=theta, X=dose_data, y=resp_data,
+    >>>                           lower=0, upper=8)
 
+    The 2PL model and others are harder to optimize compared to other estimating equations. Namely, the optimizer is
+    not aware of implicit bounds on the parameters. To reduce non-convergence issues, we can give the root-finder good
+    starting values.
 
-    Calling the M-estimation procedure (note that `init` has 3 values now, since `X.shape[1]` is equal to 3).
+    First, the ED50 should be between the lower and upper limits. Second, the sign for the steepness depends on whether
+    the response declines (positive) or the response increases (negative). Finally, some solvers may be better suited
+    to the problem, so try a few different options.
 
-    >>> mestimation = MEstimator(stacked_equations=psi, init=[0., 0., 0.,])
-    >>> mestimation.estimate()
+    Here, we use some general starting values that should perform well in many cases. For ED50, give the mid-point
+    between the maximum response and the minimum response. The initial value for steepness is more difficult. Ideally,
+    we would give a starting value of zero, but that will fail in this example. Giving a small positive starting value
+    works in this example. Finally, we use the ``lm`` solver.
 
-    Inspecting the parameter estimates and the variance
+    Note
+    ----
+    To summarize the recommendations, be sure to examine your data (e.g., scatterplot). This will help to determine the
+    initial starting values for the root-finding procedure. Otherwise, you may come across a convergence error.
 
-    >>> mestimation.theta
-    >>> mestimation.variance
+    >>> estr = MEstimator(psi, init=[(np.max(resp_data)+np.min(resp_data)) / 2,
+    >>>                              (np.max(resp_data)+np.min(resp_data)) / 2])
+    >>> estr.estimate(solver='lm')
 
-    References
-    ----------
-    Boos DD, & Stefanski LA. (2013). M-estimation (estimating equations). In Essential Statistical Inference
-    (pp. 297-337). Springer, New York, NY.
-    """
-    X = np.asarray(X)                       # Convert to NumPy array
-    y = np.asarray(y)[:, None]              # Convert to NumPy array and ensure correct shape for matrix algebra
-    beta = np.asarray(theta)[:, None]       # Convert to NumPy array and ensure correct shape for matrix algebra
+    Inspecting the parameter estimates, variance, and confidence intervals
 
-    # Allowing for a weighted linear model
-    if weights is None:                     # If weights is unspecified
-        w = np.ones(X.shape[0])                 # ... assign weight of 1 to all observations
-    else:                                   # Otherwise
-        w = np.asarray(weights)                 # ... set weights as input vector
+    >>> estr.theta
+    >>> estr.variance
+    >>> estr.confidence_intervals()
 
-    # Output b-by-n matrix
-    return w*((y -                            # Speedy matrix algebra for regression
-               np.dot(X, beta))               # ... linear regression requires no transfomrations
-               * X).T                         # ... multiply by coefficient and transpose for correct orientation
+    Inspecting the parameter estimates
 
+    >>> estr.theta[0]    # ED(50)
+    >>> estr.theta[1]    # steepness
 
-def ee_logistic_regression(theta, X, y, weights=None):
-    r"""Default stacked estimating equation for logistic regression. The estimating equation is
+    References
+    ----------
+    Ritz C, Baty F, Streibig JC, & Gerhard D. (2015). Dose-response analysis using R. *PloS One*, 10(12), e0146021.
 
-    .. math::
+    An H, Justin TL, Aubrey GB, Marron JS, & Dittmer DP. (2019). dr4pl: A Stable Convergence Algorithm for the 4
+    Parameter Logistic Model. *R J.*, 11(2), 171.
 
-        \sum_i^n \psi(Y_i, X_i, \theta) = \sum_i^n (Y_i - expit(X_i^T \theta)) X_i = 0
+    Inderjit, Streibig JC, & Olofsdotter M. (2002). Joint action of phenolic acid mixtures and its significance in
+    allelopathy research. *Physiologia Plantarum*, 114(3), 422-428.
+    """
+    # Creating rho to cut down on typing
+    rho = (X / theta[0])**theta[1]
 
-    where expit, or the inverse logit is
+    # Generalized 3PL model function for y-hat
+    fx = lower + (upper - lower) / (1 + rho)
 
-    .. math::
+    # Using a special implementatin of natural log here
+    nested_log = np.log(X / theta[0],             # ... to avoid dose=0 issues only take log
+                        where=0 < X)              # ... where dose>0 (otherwise puts zero in place)
 
-        expit(u) = 1 / (1 + exp(u))
+    # Calculate the derivatives for the gradient
+    deriv = np.array(((upper-lower)*theta[1]/theta[0]*rho/(1+rho)**2,     # Gradient for steepness
+                      (upper-lower) * nested_log * rho / (1+rho)**2), )   # Gradient for ED50
 
-    Here, theta is a 1-by-b array, where b is the distinct covariates included as part of X. For example, if X is a
-    3-by-n matrix, then theta will be a 1-by-3 array. The code is general to allow for an arbitrary number of X's (as
-    long as there is enough support in the data).
+    # Compute gradient and return for each i
+    return -2*(y-fx)*deriv
 
-    Note
-    ----
-    All provided estimating equations are meant to be wrapped inside a user-specified function. Throughtout, these
-    user-defined functions are defined as `psi`.
 
+def ee_effective_dose_delta(theta, y, delta, steepness, ed50, lower, upper):
+    r"""Default stacked estimating equation to pair with the 4 parameter logistic model for estimation of the
+    :math:`delta` effective dose. The estimating equation is
 
-    Here, theta corresponds to the coefficients in a logistic regression model, and therefore are the log-odds.
+    .. math::
+
+        \psi(Y_i, \theta) = \beta_1 + \frac{\beta_4 - \beta_1}{1 + (\theta / \beta_2)^{\beta_3}} - \beta_4(1-\delta)
+        - \beta_1 \delta
+
+    where theta is the :math:`ED(\delta)`, and the beta values are from a 4PL model (1: lower limit, 2: steepness,
+    3: ED(50), 4: upper limit). When lower or upper limits are place, the corresponding beta's are replaced by
+    constants. For proper uncertainty estimation, this estimating equation should be stacked together with the
+    correspond PL model.
 
     Note
     ----
-    For complex regression problems, the optimizer behind the scenes is not particularly robust (unlike functions
-    specializing in solely logistic regression). Therefore, optimization of logistic regression via a separate
-    functionality can be done then those estimated parameters are fed forward as the initial values (which should
-    result in a more stable optimization).
-
+    This estimating equation is meant to be paired with the estimating equations for either the 4PL, 3PL, or 2PL models.
 
     Parameters
     ----------
-    theta : 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 accomplished generally by `[0, ] * X.shape[1]`.
-    X : vector
-        2-dimensional vector of n observed values for b variables. No missing data should be included (missing data
-        may cause unexpected behavior).
-    y : vector
-        1-dimensional vector of n observed values. The Y values should all be 0 or 1. No missing data should be
-        included (missing data may cause unexpected behavior).
-    weights : vector, None, optional
-        1-dimensional vector of n weights. No missing weights should be included. Default is None, which assigns a
-        weight of 1 to all observations.
+    theta : int, float
+        Theta value corresponding to the ED(alpha).
+    y : ndarray, list, vector
+        1-dimensional vector of n response values, used to construct correct shape for output.
+    delta : float
+        The effective dose level of interest, ED(alpha).
+    steepness : float
+        Estimated parameter for the steepness from the PL.
+    ed50 : float
+        Estimated parameter for the ED50, or ED(alpha=50), from the PL.
+    lower : int, float
+        Estimated parameter or pre-specified constant for the lower limit. This should be a pre-specified constant for
+        both the 3PL and 2PL.
+    upper : int, float
+        Estimated parameter or pre-specified constant for the lower limit. This should be a pre-specified constant for
+        the 2PL.
 
     Returns
     -------
     array :
-        Returns a b-by-n NumPy array evaluated for the input theta and y
+        Returns a 1-by-n NumPy array evaluated for the input theta
 
     Examples
     --------
-    Construction of a estimating equation(s) with `ee_logistic_regression` should be done similar to the following
+    Construction of a estimating equations for ED25 with ``ee_3p_logistic`` 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_logistic_regression
+    >>> from delicatessen.data import load_inderjit
+    >>> from delicatessen.estimating_equations import ee_2p_logistic, ee_effective_dose_delta
 
-    Some generic data to estimate a logistic regresion model
-
-    >>> n = 500
-    >>> data = pd.DataFrame()
-    >>> data['X'] = np.random.normal(size=n)
-    >>> data['Z'] = np.random.normal(size=n)
-    >>> data['Y'] = np.random.binomial(n=1, p=logistic.cdf(0.5 + 2*data['X'] - 1*data['Z']), size=n)
-    >>> data['C'] = 1
+    For demonstration, we use dose-response data from Inderjit et al. (2002), which can be loaded from ``delicatessen``
+    directly.
 
-    Note that `C` here is set to all 1's. This will be the intercept in the regression.
+    >>> d = load_inderjit()   # Loading array of data
+    >>> dose_data = d[:, 1]   # Dose data
+    >>> resp_data = d[:, 0]   # Response data
 
-    Defining psi, or the stacked estimating equations
+    Since there is a natural lower-bound of 0 for root growth, we set ``lower=0``. While a natural upper bound does not
+    exist for this example, we set ``upper=8`` for illustrative purposes. Defining psi, or the stacked estimating
+    equations
 
     >>> def psi(theta):
-    >>>         return ee_logistic_regression(theta=theta, X=data[['C', 'X', 'Z']], y=data['Y'])
+    >>>     pl_model = ee_3p_logistic(theta=theta, X=dose_data, y=resp_data,
+    >>>                               lower=0)
+    >>>     ed_25 = ee_effective_dose_delta(theta[3], y=resp_data, delta=0.20,
+    >>>                                     steepness=theta[0], ed50=theta[1],
+    >>>                                     lower=0, upper=theta[2])
+    >>>     # Returning stacked estimating equations
+    >>>     return np.vstack((pl_model,
+    >>>                       ed_25,))
 
+    Notice that the estimating equations are stacked in the order of the parameters in ``theta`` (the first 3 belong to
+    3PL and the last belong to ED(25)).
 
-    Calling the M-estimation procedure (note that `init` has 3 values now, since `X.shape[1]` is equal to 3).
+    >>> estr = MEstimator(psi, init=[(np.max(resp_data)+np.min(resp_data)) / 2,
+    >>>                              (np.max(resp_data)+np.min(resp_data)) / 2,
+    >>>                              np.max(resp_data),
+    >>>                              (np.max(resp_data)+np.min(resp_data)) / 2])
+    >>> estr.estimate(solver='lm')
 
-    >>> mestimation = MEstimator(stacked_equations=psi, init=[0., 0., 0.,])
-    >>> mestimation.estimate()
+    Inspecting the parameter estimates, variance, and confidence intervals
 
-    Inspecting the parameter estimates and the variance
+    >>> estr.theta
+    >>> estr.variance
+    >>> estr.confidence_intervals()
 
-    >>> mestimation.theta
-    >>> mestimation.variance
+    Inspecting the parameter estimates
+
+    >>> estr.theta[0]    # ED(50)
+    >>> estr.theta[1]    # steepness
+    >>> estr.theta[2]    # upper limit
+    >>> estr.theta[3]    # ED(25)
 
     References
     ----------
-    Boos DD, & Stefanski LA. (2013). M-estimation (estimating equations). In Essential Statistical Inference
-    (pp. 297-337). Springer, New York, NY.
+    Ritz C, Baty F, Streibig JC, & Gerhard D. (2015). Dose-response analysis using R. *PloS One*, 10(12), e0146021.
+
+    An H, Justin TL, Aubrey GB, Marron JS, & Dittmer DP. (2019). dr4pl: A Stable Convergence Algorithm for the 4
+    Parameter Logistic Model. *R J.*, 11(2), 171.
+
+    Inderjit, Streibig JC, & Olofsdotter M. (2002). Joint action of phenolic acid mixtures and its significance in
+    allelopathy research. *Physiologia Plantarum*, 114(3), 422-428.
     """
-    X = np.asarray(X)                      # Convert to NumPy array
-    y = np.asarray(y)[:, None]             # Convert to NumPy array and ensure correct shape for matrix algebra
-    beta = np.asarray(theta)[:, None]      # Convert to NumPy array and ensure correct shape for matrix algebra
+    # Creating rho to cut down on typing
+    rho = (theta / steepness)**ed50            # Theta is the corresponds ED(alpha) value
 
-    # Allowing for a weighted linear model
-    if weights is None:                     # If weights is unspecified
-        w = np.ones(X.shape[0])                 # ... assign weight of 1 to all observations
-    else:                                   # Otherwise
-        w = np.asarray(weights)                 # ... set weights as input vector
+    # Calculating the predicted value for f(x,\theta), or y-hat
+    fx = lower + (upper - lower) / (1 + rho)
 
-    # Output b-by-n matrix
-    return w*((y -                                # Speedy matrix algebra for regression
-               inverse_logit(np.dot(X, beta)))    # ... inverse-logit transformation of predictions
-              * X).T                              # ... multiply by coefficient and transpose for correct orientation
+    # Subtracting off (Upper*(1-delta) + Lower*delta) since theta should result in zeroing of quantity
+    ed_delta = fx - upper*(1-delta) - lower*delta
+
+    # Returning constructed 1-by-ndarray for stacked estimating equations
+    return np.ones(np.asarray(y).shape[0])*ed_delta
 
 
 #################################################################
 # Causal Inference (ATE) Estimating Equations
 
 
-def ee_gformula(theta, X, y, treat_index, force_continuous=False):
-    r"""Default stacked estimating equation for the parametric g-formula in the time-fixed setting. The parameter(s) of
-    interest is the average treatment effect, with potential interest in the underlying risk or mean functions. For
-    continuous Y, the linear regression estimating equation is
+def ee_gformula(theta, y, X, X1, X0=None, force_continuous=False):
+    r"""Default stacked estimating equation for parametric g-computation in the time-fixed setting. The parameter of
+    interest can either be the mean under a single interventions or plans on an action, or the mean difference between
+    two interventions or plans on an action. This is accomplished by providing the estimating equation the observed
+    data (``X``, ``y``), and the same data under the actions (``X1`` and optionally ``X0``).
+
+    For continuous Y, the linear regression estimating equation is
 
     .. math::
 
@@ -454,75 +1344,88 @@ def ee_gformula(theta, X, y, treat_index, force_continuous=False):
 
     .. math::
 
-        \sum_i^n \psi_m(Y_i, X_i, \theta) = \sum_i^n (Y_i - expit(X_i^T \theta)) X_i = 0
+        \sum_i^n \psi_m(Y_i, X_i, \beta) = \sum_i^n (Y_i - expit(X_i^T \beta)) X_i = 0
 
     By default, `ee_gformula` detects whether `y` is all binary (zero or one), and applies logistic regression if that
-    happens. See the parameters for more details.
+    is evaluated to be true. See the parameters for further details.
 
-    For the implementation of the g-formula, stacked estimating equations are also used for the risk / mean had
-    everyone been given treatment=1, the risk / mean had everyone been given treatment=0, and the risk / mean
-    difference between those two risks. Respectively, those estimating equations look like
+    There are two variations on the parameter of interest. The first could be the mean under a plan, where the plan sets
+    the values of action :math:`A` (e.g., exposure, treatment, vaccination, etc.). The estimating equation for this
+    causal mean is
 
     .. math::
 
         \sum_i^n \psi_1(Y_i, X_i, \theta_1) = \sum_i^n g(\hat{Y}_i) - \theta_1 = 0
 
-    .. math::
+    Here, the function :math:`g(.)` is a generic function. If linear regression was used, :math:`g(.)` is the identity
+    function. If logistic regression was used, :math:`g(.)` is the expit or inverse-logit function.
 
-        \sum_i^n \psi_0(Y_i, X_i, \theta_2) = \sum_i^n g(\hat{Y}_i) - \theta_2 = 0
+    Note
+    ----
+    This variation includes :math:`1+b` parameters, where the first parameter is the causal mean, and the remainder are
+    the parameters for the regression model.
+
+    The alternative parameter of interest could be the mean difference between two plans. A common example of this would
+    be the average causal effect, where the plans are all-action-one versus all-action-zero. Therefore, the estimating
+    equations consist of the following three equations
 
     .. math::
 
         \sum_i^n \psi_0(Y_i, X_i, \theta_0) = \sum_i^n (\theta_1 - \theta_2) - \theta_0 = 0
 
-    Here, the function g() is a generic function. If linear regression was used, g() is the identity function. If
-    logistic regression was used, g() is the expit or inverse-logit function.
+        \sum_i^n \psi_1(Y_i, X_i, \theta_1) = \sum_i^n g(\hat{Y}_i) - \theta_1 = 0
+
+        \sum_i^n \psi_0(Y_i, X_i, \theta_2) = \sum_i^n g(\hat{Y}_i) - \theta_2 = 0
 
-    Due to these 3 extra values, the length of the theta vector is b+3, where b is the number of parameters in the
-    regression model.
 
     Note
     ----
-    All provided estimating equations are meant to be wrapped inside a user-specified function. Throughtout, these
-    user-defined functions are defined as `psi`.
-
+    This variation includes :math:`3+b` parameters, where the first parameter is the causal mean difference, the second
+    is the causal mean under plan 1, the third is the causal mean under plan 0, and the remainder are the parameters
+    for the regression model.
 
-    Here, theta corresponds to a variety of different quantities. The *first* value in theta vector is the risk / mean
-    difference (or average treatment effect), the *second* is the risk / mean had everyone been given treatment=0, the
-    *third* is the risk / mean had everyone been given treatment=1. The remainder of the parameters correspond to the
-    regression model coefficients, in the order input.
+    The parameter of interest is designated by the user via whether the optional argument ``X0`` is left as ``None``
+    (which estimates the causal mean) or is given an array (which estimates the causal mean difference and the
+    corresponding causal means).
 
     Note
     ----
-    For complex regression problems, the optimizer behind the scenes is not particularly robust (unlike functions
-    specializing in solely regression models). Therefore, optimization of the regression model via a separate
-    functionality can be done then those estimated parameters are fed forward as the initial values (which should
-    result in a more stable optimization).
+    All provided estimating equations are meant to be wrapped inside a user-specified function. Throughtout, these
+    user-defined functions are defined as ``psi``.
 
+    See the examples below for how action plans are specified.
 
     Parameters
     ----------
-    theta : array, list
+    theta : ndarray, list, vector
         Array of parameters to estimate. For the Cox model, corresponds to the log hazard ratios
-    X : vector
-        2-dimensional vector of n observed values for b variables. No missing data should be included (missing data
-        may cause unexpected behavior).
-    y : vector
+    y : ndarray, list, vector
         1-dimensional vector of n observed values. The Y values should all be 0 or 1. No missing data should be
         included (missing data may cause unexpected behavior).
-    treat_index : int
-        Column index for the treatment vector.
+    X : ndarray, list, vector
+        2-dimensional vector of n observed values for b variables. No missing data should be included (missing data
+        may cause unexpected behavior).
+    X1 : ndarray, list, vector
+        2-dimensional vector of n observed values for b variables under the action plan. If the action is indicated by
+        ``A``, then ``X1`` will take the original data ``X`` and update the values of ``A`` to follow the deterministic
+        plan. No missing data should be included (missing data may cause unexpected behavior).
+    X0 : ndarray, list, vector, None, optional
+        2-dimensional vector of n observed values for b variables under the action plan. This second argument is
+        optional and should be specified if a causal mean difference between two action plans is of interest. If the
+        action is indicated by ``A``, then ``X0`` will take the original data ``X`` and update the values of ``A`` to
+        follow the deterministic reference plan. No missing data should be included (missing data may cause unexpected
+        behavior).
     force_continuous : bool, optional
         Option to force the use of linear regression despite detection of a binary variable.
 
     Returns
     -------
     array :
-        Returns a (3+b)-by-n NumPy array evaluated for the input theta and y
+        Returns a (1+b)-by-n NumPy array if ``X0=None``, or returns a (3+b)-by-n NumPy array if ``X0!=None``
 
     Examples
     --------
-    Construction of a estimating equation(s) with `ee_gformula` should be done similar to the following
+    Construction of a estimating equation(s) with ``ee_gformula`` should be done similar to the following
 
     >>> import numpy as np
     >>> import pandas as pd
@@ -540,153 +1443,202 @@ def ee_gformula(theta, X, y, treat_index, force_continuous=False):
     >>> d['Y'] = (1-d['A'])*d['Ya0'] + d['A']*d['Ya1']
     >>> d['C'] = 1
 
-    Defining psi, or the stacked estimating equations. Note that 'A' is the treatment of interest, so `treat_index` is
-    set to 1 (compared to input `X`).
+    In the first example, we will estimate the causal mean had everyone been set to ``A=1``. Therefore, the optional
+    argument ``X0`` is left as ``None``. Before creating the estimating equation, we need to do some data prep. First,
+    we will create an interaction term between ``A`` and ``W`` in the original data. Then we will generate a copy of
+    the data and update the values of ``A`` to be all ``1``.
+
+    >>> d['AW'] = d['A']*d['W']
+    >>> d1 = d.copy()
+    >>> d1['A'] = 1
+    >>> d1['AW'] = d1['A']*d1['W']
+
+    Having setup our data, we can now define the psi function.
 
     >>> def psi(theta):
-    >>>     return ee_gformula(theta, X=d[['C', 'A', 'W']], y=d['Y'], treat_index=1)
+    >>>     return ee_gformula(theta,
+    >>>                        y=d['Y'],
+    >>>                        X=d[['C', 'A', 'W', 'AW']],
+    >>>                        X1=d1[['C', 'A', 'W', 'AW']])
+
+    Notice that ``y`` corresponds to the observed outcomes, ``X`` corresponds to the observed covariate data, and ``X1``
+    corresponds to the covariate data *under the action plan*.
+
+    Now we can call the M-Estimation procedure. Since we are estimating the causal mean, and the regression parameters,
+    the length of the initial values needs to correspond with this. Our linear regression model consists of 4
+    coefficients, so we need 1+4=5 initial values. When the outcome is binary (like it is in this example), we can be
+    nice to the optimizer and give it a starting value of 0.5 for the causal mean (since 0.5 is in the middle of that
+    distribution). Below is the call to ``MEstimator``
+
+    >>> estr = MEstimator(psi, init=[0.5, 0., 0., 0., 0.])
+    >>> estr.estimate(solver='lm')
+
+    Inspecting the parameter estimates, variance, and 95% confidence intervals
 
-    Calling the M-estimation procedure. Since `X` is 3-by-n here and g-formula has 3 additional parameters, the initial
-    values should be of length 3+3=6. In general, it will be best to start with [0., 0.5, 0.5, ...] as the initials for
-    the risk parameters. This will start the initial at the exact middle value for each of the first 3 parameters. For
-    the regression coefficients, those can be set as zero, or if there is difficulty in simultaneous optimization,
-    coefficient estimates from outside `MEstimator` can be provided as inputs.
+    >>> estr.theta
+    >>> estr.variance
+    >>> estr.confidence_intervals()
 
-    >>> mestimation = MEstimator(stacked_equations=psi, init=[0., 0.5, 0.5, 0., 0., 0.])
-    >>> mestimation.estimate()
+    More specifically, the causal mean is
 
-    Inspecting the parameter estimates and the variance
+    >>> estr.theta[0]
 
-    >>> mestimation.theta
-    >>> mestimation.variance
+    Continuing from the previous example, let's say we wanted to estimate the average causal effect. Therefore, we want
+    to contrast two plans (all ``A=1`` versus all ``A=0``). As before, we need to create the reference data for ``X0``
 
-    More specifically, the average treatment effect and its variance are
+    >>> d0 = d.copy()
+    >>> d0['A'] = 0
+    >>> d0['AW'] = d0['A']*d0['W']
 
-    >>> mestimation.theta[0]
-    >>> mestimation.variance[0, 0]
+    Having setup our data, we can now define the psi function.
 
-    The risk / mean had all been given treatment=1
+    >>> def psi(theta):
+    >>>     return ee_gformula(theta,
+    >>>                        y=d['Y'],
+    >>>                        X=d[['C', 'A', 'W', 'AW']],
+    >>>                        X1=d1[['C', 'A', 'W', 'AW']],
+    >>>                        X0=d0[['C', 'A', 'W', 'AW']], )
+
+    Notice that ``y`` corresponds to the observed outcomes, ``X`` corresponds to the observed covariate data, ``X1``
+    corresponds to the covariate data under ``A=1``, and ``X0`` corresponds to the covariate data under ``A=0``. Here,
+    we need 3+4=7 starting values, since there are two additional parameters from the previous example. For the
+    difference, a starting value of 0 is generally a good choice. Since ``Y`` is binary, we again provide 0.5 as
+    starting values for the causal means
 
-    >>> mestimation.theta[1]
-    >>> mestimation.variance[1, 1]
+    >>> estr = MEstimator(psi, init=[0., 0.5, 0.5, 0., 0., 0., 0.])
+    >>> estr.estimate(solver='lm')
 
-    The risk / mean had all been given treatment=0
+    Inspecting the parameter estimates
 
-    >>> mestimation.theta[2]
-    >>> mestimation.variance[2, 2]
+    >>> estr.theta[0]    # causal mean difference of 1 versus 0
+    >>> estr.theta[1]    # causal mean under X1
+    >>> estr.theta[2]    # causal mean under X0
+    >>> estr.theta[3:]   # logistic regression coefficients
 
     References
     ----------
     Snowden JM, Rose S, & Mortimer KM. (2011). Implementation of G-computation on a simulated data set: demonstration
-    of a causal inference technique. American Journal of Epidemiology, 173(7), 731-738.
+    of a causal inference technique. *American Journal of Epidemiology*, 173(7), 731-738.
+
+    Hernán MA, & Robins JM. (2006). Estimating causal effects from epidemiological data.
+    *Journal of Epidemiology & Community Health*, 60(7), 578-586.
     """
     # Ensuring correct typing
-    X = np.asarray(X)                            # Convert to NumPy array
-    y = np.asarray(y)                            # Convert to NumPy array
-    beta = theta[3:]                       # Extracting out theta's for the regression model
+    X = np.asarray(X)                # Convert to NumPy array
+    y = np.asarray(y)                # Convert to NumPy array
+    X1 = np.asarray(X1)              # Convert to NumPy array
+
+    # Error checking for misaligned shapes
+    if X.shape != X1.shape:
+        raise ValueError("The dimensions of X and X1 must be the same.")
+
+    # Processing data depending on if two plans were specified
+    if X0 is None:                   # If no reference was specified
+        mu1 = theta[0]                  # ... only a single mean
+        beta = theta[1:]                # ... immediately followed by the regression parameters
+    else:                            # Otherwise difference and both plans are to be returned
+        X0 = np.asarray(X0)             # ... reference data to NumPy array
+        if X.shape != X0.shape:         # ... error checking for misaligned shapes
+            raise ValueError("The dimensions of X and X0 must be the same.")
+        mud = theta[0]                  # ... first parameter is mean difference
+        mu1 = theta[1]                  # ... second parameter is mean under X1
+        mu0 = theta[2]                  # ... third parameter is mean under X0
+        beta = theta[3:]                # ... remainder are for the regression model
 
     # Checking outcome variable type
     if np.isin(y, [0, 1]).all() and not force_continuous:
-        continuous = False
+        regression = ee_logistic_regression         # Use a logistic regression model
+        transform = inverse_logit                   # ... and need to inverse-logit transformation
     else:
-        continuous = True
-
-    if continuous:
-        # Linear regression parameters
-        preds_reg = ee_linear_regression(theta=beta,
-                                         X=X,
-                                         y=y)
-        # Calculating Y(a=0)
-        X[:, treat_index] = 0
-        ya0 = np.dot(X, beta) - theta[2]
-        # Calculating Y(a=1)
-        X[:, treat_index] = 1
-        ya1 = np.dot(X, beta) - theta[1]
-    else:
-        # Logistic regression parameters
-        preds_reg = ee_logistic_regression(theta=beta,
-                                           X=X,
-                                           y=y)
-        # Calculating Y(a=0)
-        X[:, treat_index] = 0
-        ya0 = inverse_logit(np.dot(X, beta)) - theta[2]
-        # Calculating Y(a=1)
-        X[:, treat_index] = 1
-        ya1 = inverse_logit(np.dot(X, beta)) - theta[1]
-
-    # Calculating Y(a=1) - Y(a=0)
-    ate = np.ones(y.shape[0]) * (theta[1] - theta[2]) - theta[0]
-
-    # Output (3+b)-by-n stacked array
-    return np.vstack((ate,            # theta[0] is for the ATE
-                      ya1[None, :],   # theta[1] is for R1
-                      ya0[None, :],   # theta[2] is for R0
-                      preds_reg))     # theta[3:] is for the regression coefficients
+        regression = ee_linear_regression           # Use a linear regression model
+        transform = identity                        # ... and need to apply the identity (no) transformation
+
+    # Estimating regression parameters
+    preds_reg = regression(theta=beta,              # beta coefficients
+                           X=X, y=y)                # along with observed X and observed y
+
+    # Calculating mean under X1
+    ya1 = transform(np.dot(X1, beta)) - mu1         # mean under X1
+
+    if X0 is None:                                  # if no X0, then nothing left to do
+        # Output (1+b)-by-n stacked array
+        return np.vstack((ya1[None, :],     # theta[0] is the mean under X1
+                          preds_reg))       # theta[1:] is the regression coefficients
+    else:                                           # if X0, then need to predict mean under X0 and difference
+        # Calculating mean under X0
+        ya0 = transform(np.dot(X0, beta)) - mu0
+        # Calculating mean difference between X1 and X0
+        ace = np.ones(y.shape[0])*(mu1 - mu0) - mud
+        # Output (3+b)-by-n stacked array
+        return np.vstack((ace,            # theta[0] is the mean difference between X1 and X0
+                          ya1[None, :],   # theta[1] is the mean under X1
+                          ya0[None, :],   # theta[2] is the mean under X0
+                          preds_reg))     # theta[3:] is for the regression coefficients
+
+
+def ee_ipw(theta, y, A, W, truncate=None):
+    r"""Default stacked estimating equation for inverse probability weighting in the time-fixed setting. The
+    parameter of interest is the average causal effect. For estimation of the weights (or propensity scores), a
+    logistic model is used.
 
+    Note
+    ----
+    Unlike ``ee_gformula``, ``ee_ipw`` only provides the average causal effect (and the causal means for ``A=1`` and
+    ``A=0``). In other words, the implementation of IPW does not support generic action plans off-the-shelf,
+    unlike ``ee_gformula``.
 
-def ee_ipw(theta, X, y, treat_index):
-    r"""Default stacked estimating equation for inverse probability weighting in the time-fixed setting. The
-    parameter(s) of interest is the average treatment effect, with potential interest in the underlying risk or mean
-    functions. For estimation of the weights (or propensity scores), a logistic model is used. Therefore, the first
-    estimating equation is
+    The first estimating equation for the logistic regression model is
 
     .. math::
 
-        \sum_i^n \psi_g(A_i, W_i, \theta) = \sum_i^n (A_i - expit(W_i^T \theta)) W_i = 0
+        \sum_i^n \psi_g(A_i, W_i, \alpha) = \sum_i^n (A_i - expit(W_i^T \alpha)) W_i = 0
 
-    where A is the treatment and W is the set of confounders. Both of these are processed from the input `X` and the
-    specified `treat_index`.
+    where A is the treatment and W is the set of confounders.
 
-    For the implementation of the inverse probability weighting estimator, stacked estimating equations are also used
-    for the risk / mean had everyone been given treatment=1, the risk / mean had everyone been given treatment=0, and
-    the risk / mean difference between those two risks. Respectively, those estimating equations look like
+    For the implementation of the inverse probability weighting estimator, stacked estimating equations are used
+    for the mean had everyone been set to ``A=1``, the mean had everyone been set to ``A=0``, and the mean difference
+    between the two causal means. The estimating equations are
 
     .. math::
 
-        \sum_i^n \psi_1(Y_i, A_i, \pi_i, \theta_1) = \sum_i^n \frac{A_i \times Y_i}{\pi_i} - \theta_1 = 0
+        \sum_i^n \psi_d(Y_i, A_i, \pi_i, \theta_0) = \sum_i^n (\theta_1 - \theta_2) - \theta_0 = 0
 
-    .. math::
+        \sum_i^n \psi_1(Y_i, A_i, \pi_i, \theta_1) = \sum_i^n \frac{A_i \times Y_i}{\pi_i} - \theta_1 = 0
 
         \sum_i^n \psi_0(Y_i, A_i, \pi_i, \theta_2) = \sum_i^n \frac{(1-A_i) \times Y_i}{1-\pi_i} - \theta_2 = 0
 
-    .. math::
-
-        \sum_i^n \psi_d(Y_i, A_i, \pi_i, \theta_0) = \sum_i^n (\theta_1 - \theta_2) - \theta_0 = 0
 
-    Due to these 3 extra values, the length of the theta vector is b+3, where b is the number of parameters in the
+    Due to these 3 extra values, the length of the theta vector is 3+b, where b is the number of parameters in the
     regression model.
 
     Note
     ----
     All provided estimating equations are meant to be wrapped inside a user-specified function. Throughtout, these
-    user-defined functions are defined as `psi`.
-
-
-    Here, theta corresponds to a variety of different quantities. The *first* value in theta vector is the risk / mean
-    difference (or average treatment effect), the *second* is the risk / mean had everyone been given treatment=0, the
-    *third* is the risk / mean had everyone been given treatment=1. The remainder of the parameters correspond to the
-    logistic regression model coefficients, in the order input.
-
-    Note
-    ----
-    For complex regression problems, the optimizer behind the scenes is not particularly robust (unlike functions
-    specializing in solely regression models). Therefore, optimization of the regression model via a separate
-    functionality can be done then those estimated parameters are fed forward as the initial values (which should
-    result in a more stable optimization).
+    user-defined functions are defined as ``psi``.
 
+    Here, theta corresponds to a variety of different quantities. The *first* value in theta vector is the mean
+    difference (or average causal effect), the *second* is the mean had everyone been set to ``A=1``, the *third* is the
+    mean had everyone been set to ``A=0``. The remainder of the parameters correspond to the logistic regression model
+    coefficients.
 
     Parameters
     ----------
-    theta : array, list
+    theta : ndarray, list, vector
         Array of parameters to estimate. For the Cox model, corresponds to the log hazard ratios
-    X : vector
-        2-dimensional vector of n observed values for b variables. No missing data should be included (missing data
-        may cause unexpected behavior).
-    y : vector
-        1-dimensional vector of n observed values. The Y values should all be 0 or 1. No missing data should be
+    y : ndarray, list, vector
+        1-dimensional vector of n observed values. No missing data should be included (missing data may cause
+        unexpected behavior).
+    A : ndarray, list, vector
+        1-dimensional vector of n observed values. The A values should all be 0 or 1. No missing data should be
         included (missing data may cause unexpected behavior).
-    treat_index : int
-        Column index for the treatment vector.
+    W : ndarray, list, vector
+        2-dimensional vector of n observed values for b variables to model the probability of ``A`` with. No missing
+        data should be included (missing data may cause unexpected behavior).
+    truncate : None, list, set, optional
+        Bounds to truncate the estimated probabilities of ``A`` at. For example, estimated probabilities above 0.99 or
+        below 0.01 can be set to 0.99 or 0.01, respectively. This is done by specifying ``truncate=(0.01, 0.99)``. Note
+        this step is done via ``numpy.clip(.., a_min=truncate[0], a_max=truncate[1])``, so order is important. Default
+        is None, which applies to no truncation.
 
     Returns
     -------
@@ -695,7 +1647,7 @@ def ee_ipw(theta, X, y, treat_index):
 
     Examples
     --------
-    Construction of a estimating equation(s) with `ee_ipw` should be done similar to the following
+    Construction of a estimating equation(s) with ``ee_ipw`` should be done similar to the following
 
     >>> import numpy as np
     >>> import pandas as pd
@@ -713,54 +1665,49 @@ def ee_ipw(theta, X, y, treat_index):
     >>> d['Y'] = (1-d['A'])*d['Ya0'] + d['A']*d['Ya1']
     >>> d['C'] = 1
 
-    Defining psi, or the stacked estimating equations. Note that 'A' is the treatment of interest, so `treat_index` is
-    set to 1 (compared to input `X`).
+    Defining psi, or the stacked estimating equations. Note that 'A' is the action.
 
     >>> def psi(theta):
-    >>>     return ee_ipw(theta, X=d[['C', 'A', 'W']], y=d['Y'], treat_index=1)
+    >>>     return ee_ipw(theta, y=d['Y'], A=d['A'],
+    >>>                   W=d[['C', 'W']])
 
-    Calling the M-estimation procedure. Since `X` is 3-by-n here and g-formula has 3 additional parameters, the initial
-    values should be of length 3+3=6. In general, it will be best to start with [0., 0.5, 0.5, ...] as the initials for
-    the risk parameters. This will start the initial at the exact middle value for each of the first 3 parameters. For
-    the regression coefficients, those can be set as zero, or if there is difficulty in simultaneous optimization,
-    coefficient estimates from outside `MEstimator` can be provided as inputs.
+    Calling the M-estimation procedure. Since `X` is 2-by-n here and IPW has 3 additional parameters, the initial
+    values should be of length 3+2=5. In general, it will be best to start with [0., 0.5, 0.5, ...] as the initials when
+    ``Y`` is binary. Otherwise, starting with all 0. as initials is reasonable.
 
-    >>> mestimation = MEstimator(stacked_equations=psi, init=[0., 0.5, 0.5, 0., 0., 0.])
-    >>> mestimation.estimate()
+    >>> estr = MEstimator(stacked_equations=psi, init=[0., 0.5, 0.5, 0., 0.])
+    >>> estr.estimate(solver='lm')
 
-    Inspecting the parameter estimates and the variance
+    Inspecting the parameter estimates, variance, and 95% confidence intervals
 
-    >>> mestimation.theta
-    >>> mestimation.variance
+    >>> estr.theta
+    >>> estr.variance
+    >>> estr.confidence_intervals()
 
-    More specifically, the average treatment effect and its variance are
+    More specifically, the corresponding parameters are
 
-    >>> mestimation.theta[0]
-    >>> mestimation.variance[0, 0]
+    >>> estr.theta[0]    # causal mean difference of 1 versus 0
+    >>> estr.theta[1]    # causal mean under A=1
+    >>> estr.theta[2]    # causal mean under A=0
+    >>> estr.theta[3:]   # logistic regression coefficients
 
-    The risk / mean had all been given treatment=1
-
-    >>> mestimation.theta[1]
-    >>> mestimation.variance[1, 1]
-
-    The risk / mean had all been given treatment=0
-
-    >>> mestimation.theta[2]
-    >>> mestimation.variance[2, 2]
+    If you want to see how truncating the probabilities works, try repeating the above code but specifying
+    ``truncate=(0.1, 0.9)`` as an optional argument in ``ee_ipw``.
 
     References
     ----------
-    ... looking for a good one ...
+    Hernán MA, & Robins JM. (2006). Estimating causal effects from epidemiological data.
+    *Journal of Epidemiology & Community Health*, 60(7), 578-586.
+
+    Cole SR, & Hernán MA. (2008). Constructing inverse probability weights for marginal structural models.
+    *American Journal of Epidemiology*, 168(6), 656-664.
     """
     # Ensuring correct typing
-    X = np.asarray(X)                            # Convert to NumPy array
+    W = np.asarray(W)                            # Convert to NumPy array
+    A = np.asarray(A)                            # Convert to NumPy array
     y = np.asarray(y)                            # Convert to NumPy array
     beta = theta[3:]                             # Extracting out theta's for the regression model
 
-    # Splitting X into A,W (treatment, covariates)
-    W = np.delete(X, treat_index, axis=1)        # Extract all-but treatment col A
-    A = X[:, treat_index]                        # Extract treatment col A
-
     # Estimating propensity score
     preds_reg = ee_logistic_regression(theta=beta,    # Using logistic regression
                                        X=W,           # Plug-in covariates for X
@@ -768,13 +1715,15 @@ def ee_ipw(theta, X, y, treat_index):
 
     # Estimating weights
     pi = inverse_logit(np.dot(W, beta))          # Getting Pr(A|W) from model
+    if truncate is not None:                     # Truncating Pr(A|W) when requested
+        if truncate[0] > truncate[1]:
+            raise ValueError("truncate values must be specified in ascending order")
+        pi = np.clip(pi, a_min=truncate[0], a_max=truncate[1])
 
     # Calculating Y(a=1)
     ya1 = (A * y) / pi - theta[1]                # i's contribution is (AY) / \pi
-
     # Calculating Y(a=0)
     ya0 = ((1-A) * y) / (1-pi) - theta[2]        # i's contribution is ((1-A)Y) / (1-\pi)
-
     # Calculating Y(a=1) - Y(a=0)
     ate = np.ones(y.shape[0]) * (theta[1] - theta[2]) - theta[0]
 
@@ -785,99 +1734,115 @@ def ee_ipw(theta, X, y, treat_index):
                       preds_reg))      # theta[3:] is for the regression coefficients
 
 
-def ee_aipw(theta, X, y, treat_index, force_continuous=False):
-    r"""Default stacked estimating equation for augmented inverse probability weighting in the time-fixed setting. The
-    parameter(s) of interest is the average treatment effect, with potential interest in the underlying risk or mean
-    functions. For estimation of the weights (or propensity scores), a logistic model is used. Therefore, the first
-    estimating equation is
+def ee_aipw(theta, y, A, W, X, X1, X0, truncate=None, force_continuous=False):
+    r"""Default stacked estimating equation for augmented inverse probability weighting (AIPW) in the time-fixed
+    setting. The parameter of interest is the average causal effect.
+
+    Note
+    ----
+    Unlike ``ee_gformula``, ``ee_ipw`` only provides the average causal effect (and the causal means for ``A=1`` and
+    ``A=0``). In other words, the implementation of IPW does not support generic action plans off-the-shelf,
+    unlike ``ee_gformula``.
+
+    AIPW consists of two nuisance models (the propensity score model and the outcome model). For estimation of the
+    propensity scores, a logistic model is used.
 
     .. math::
 
-        \sum_i^n \psi_g(A_i, W_i, \theta) = \sum_i^n (A_i - expit(W_i^T \theta)) W_i = 0
+        \sum_i^n \psi_g(A_i, W_i, \alpha) = \sum_i^n (A_i - expit(W_i^T \alpha)) W_i = 0
 
-    where A is the treatment and W is the set of confounders. Both of these are processed from the input `X` and the
-    specified `treat_index`.
+    where ``A`` is the treatment and ``W`` is the set of confounders.
 
     Next, an outcome model is specified. For continuous Y, the linear regression estimating equation is
 
     .. math::
 
-        \sum_i^n \psi_m(Y_i, X_i, \theta) = \sum_i^n (Y_i - X_i^T \theta) X_i = 0
+        \sum_i^n \psi_m(Y_i, X_i, \beta) = \sum_i^n (Y_i - X_i^T \beta) X_i = 0
 
     and for logistic regression, the estimating equation is
 
     .. math::
 
-        \sum_i^n \psi_m(Y_i, X_i, \theta) = \sum_i^n (Y_i - expit(X_i^T \theta)) X_i = 0
+        \sum_i^n \psi_m(Y_i, X_i, \beta) = \sum_i^n (Y_i - expit(X_i^T \beta)) X_i = 0
 
     By default, `ee_aipw` detects whether `y` is all binary (zero or one), and applies logistic regression if that
-    happens. See the parameters for more details.
+    happens. See the parameters for more details. Notice that ``X`` here should consists of both ``A`` and ``W`` (with
+    possible interaction terms or other differences in functional forms from the propensity score model).
 
-    For the implementation of the augmented inverse probability weighting estimator, stacked estimating equations are
-    also used for the risk / mean had everyone been given treatment=1, the risk / mean had everyone been given
-    treatment=0, and the risk / mean difference between those two risks. Respectively, those estimating equations
+    For the implementation of the AIPW estimator, stacked estimating equations further include the mean had everyone
+    been set to ``A=1``, the mean had everyone been set to ``A=0``, and the mean difference. Those estimating equations
     look like
 
-    .. math::
-
-        \sum_i^n \psi_1(Y_i, A_i, W_i, \pi_i, \theta_1) = \sum_i^n (\frac{A_i \times Y_i}{\pi_i} - \frac{\hat{Y^1}(A_i-\pi_i}{\pi_i}) - \theta_1 = 0
-
-    .. math::
-
-        \sum_i^n \psi_0(Y_i, A_i, \pi_i, \theta_2) = \sum_i^n (\frac{(1-A_i) \times Y_i}{1-\pi_i} + \frac{\hat{Y^0}(A_i-\pi_i}{1-\pi_i})) - \theta_2 = 0
-
     .. math::
 
         \sum_i^n \psi_0(Y_i, A_i, \pi_i, \theta_0) = \sum_i^n (\theta_1 - \theta_2) - \theta_0 = 0
 
-    where Y^a is the predicted values of Y from the outcome model under treatment assignment A=a.
+        \sum_i^n \psi_1(Y_i, A_i, W_i, \pi_i, \theta_1) = \sum_i^n (\frac{A_i \times Y_i}{\pi_i} -
+        \frac{\hat{Y^1}(A_i-\pi_i}{\pi_i}) - \theta_1 = 0
 
-    Due to these 3 extra values and two nuisance models, the length of the theta vector is b+(b-1)+3, where b is the
-    number of covariats included in `X`.
-
-    Note
-    ----
-    All provided estimating equations are meant to be wrapped inside a user-specified function. Throughtout, these
-    user-defined functions are defined as `psi`.
+        \sum_i^n \psi_0(Y_i, A_i, \pi_i, \theta_2) = \sum_i^n (\frac{(1-A_i) \times Y_i}{1-\pi_i} +
+        \frac{\hat{Y^0}(A_i-\pi_i}{1-\pi_i})) - \theta_2 = 0
 
+    where :math:`Y^a` is the predicted values of :math:`Y` from the outcome model under action
+    assignment :math:`A=a`.
 
-    Here, theta corresponds to a variety of different quantities. The *first* value in theta vector is the risk / mean
-    difference (or average treatment effect), the *second* is the risk / mean had everyone been given treatment=0, the
-    *third* is the risk / mean had everyone been given treatment=1. The remainder of the parameters correspond to the
-    regression model coefficients, in the order input. The first 'chuck' of coefficients correspond to the outcome model
-    and the last 'chuck' correspond to the treatment model.
+    Due to these 3 extra values and two nuisance models, the length of the theta vector is 3+b+c, where b is the number
+    of columns in ``W``, and c is the number of columns in ``X``.
 
     Note
     ----
-    For complex regression problems, the optimizer behind the scenes is not particularly robust (unlike functions
-    specializing in solely regression models). Therefore, optimization of the regression model via a separate
-    functionality can be done then those estimated parameters are fed forward as the initial values (which should
-    result in a more stable optimization).
+    All provided estimating equations are meant to be wrapped inside a user-specified function. Throughtout, these
+    user-defined functions are defined as ``psi``.
 
+    Here, theta corresponds to a variety of different quantities. The *first* value in theta vector is mean
+    difference (or average causal effect), the *second* is the mean had everyone been given ``A=1``, the *third* is the
+    mean had everyone been given ``A=0``. The remainder of the parameters correspond to the regression model
+    coefficients, in the order input. The first 'chunk' of coefficients correspond to the propensity score model
+    and the last 'chunk' correspond to the outcome model.
 
     Parameters
     ----------
-    theta : array, list
+    theta : ndarray, list, vector
         Array of parameters to estimate. For the Cox model, corresponds to the log hazard ratios
-    X : vector
-        2-dimensional vector of n observed values for b variables. No missing data should be included (missing data
-        may cause unexpected behavior).
-    y : vector
-        1-dimensional vector of n observed values. The Y values should all be 0 or 1. No missing data should be
+    y : ndarray, list, vector
+        1-dimensional vector of n observed values. No missing data should be included (missing data may cause
+        unexpected behavior).
+    A : ndarray, list, vector
+        1-dimensional vector of n observed values. The A values should all be 0 or 1. No missing data should be
         included (missing data may cause unexpected behavior).
-    treat_index : int
-        Column index for the treatment vector.
+    W : ndarray, list, vector
+        2-dimensional vector of n observed values for b variables to model the probability of ``A`` with. No missing
+        data should be included (missing data may cause unexpected behavior).
+    X : ndarray, list, vector
+        2-dimensional vector of n observed values for c variables to model the outcome ``y``. No missing data should
+        be included (missing data may cause unexpected behavior).
+    X1 : ndarray, list, vector
+        2-dimensional vector of n observed values for b variables under the action plan. If the action is indicated by
+        ``A``, then ``X1`` will take the original data ``X`` and update the values of ``A`` to follow the deterministic
+        plan where ``A=1`` for all observations. No missing data should be included (missing data may cause unexpected
+        behavior).
+    X0 : ndarray, list, vector, None, optional
+        2-dimensional vector of n observed values for b variables under the action plan. This second argument is
+        optional and should be specified if a causal mean difference between two action plans is of interest. If the
+        action is indicated by ``A``, then ``X0`` will take the original data ``X`` and update the values of ``A`` to
+        follow the deterministic plan where ``A=0`` for all observatons. No missing data should be included (missing
+        data may cause unexpected behavior).
+    truncate : None, list, set, optional
+        Bounds to truncate the estimated probabilities of ``A`` at. For example, estimated probabilities above 0.99 or
+        below 0.01 can be set to 0.99 or 0.01, respectively. This is done by specifying ``truncate=(0.01, 0.99)``. Note
+        this step is done via ``numpy.clip(.., a_min=truncate[0], a_max=truncate[1])``, so order is important. Default
+        is None, which applies to no truncation.
     force_continuous : bool, optional
         Option to force the use of linear regression despite detection of a binary variable.
 
     Returns
     -------
     array :
-        Returns a (3+b+(b-1))-by-n NumPy array evaluated for the input theta and y
+        Returns a (3+b+c)-by-n NumPy array evaluated for the input theta and y
 
     Examples
     --------
-    Construction of a estimating equation(s) with `ee_aipw` should be done similar to the following
+    Construction of a estimating equation(s) with ``ee_aipw`` should be done similar to the following
 
     >>> import numpy as np
     >>> import pandas as pd
@@ -895,95 +1860,116 @@ def ee_aipw(theta, X, y, treat_index, force_continuous=False):
     >>> d['Y'] = (1-d['A'])*d['Ya0'] + d['A']*d['Ya1']
     >>> d['C'] = 1
 
-    Defining psi, or the stacked estimating equations. Note that 'A' is the treatment of interest, so `treat_index` is
-    set to 1 (compared to input `X`).
-
-    >>> def psi(theta):
-    >>>     return ee_aipw(theta, X=d[['C', 'A', 'W']], y=d['Y'], treat_index=1)
-
-    Calling the M-estimation procedure. Since `X` is 3-by-n here and g-formula has 3 additional parameters, the initial
-    values should be of length 3+3=6. In general, it will be best to start with [0., 0.5, 0.5, ...] as the initials for
-    the risk parameters. This will start the initial at the exact middle value for each of the first 3 parameters. For
-    the regression coefficients, those can be set as zero, or if there is difficulty in simultaneous optimization,
-    coefficient estimates from outside `MEstimator` can be provided as inputs.
+    Defining psi, or the stacked estimating equations. Note that ``A`` is the action of interest. First, we will apply
+    some necessary data processing.  We will create an interaction term between ``A`` and ``W`` in the original data.
+    Then we will generate a copy of the data and update the values of ``A=1``, and then generate another
+    copy but set ``A=0`` in that copy.
 
-    >>> mestimation = MEstimator(stacked_equations=psi, init=[0., 0.5, 0.5, 0., 0., 0., 0., 0.])
-    >>> mestimation.estimate()
+    >>> d['AW'] = d['A']*d['W']
+    >>> d1 = d.copy()   # Copy where all A=1
+    >>> d1['A'] = 1
+    >>> d1['AW'] = d1['A']*d1['W']
+    >>> d0 = d.copy()   # Copy where all A=0
+    >>> d0['A'] = 0
+    >>> d0['AW'] = d0['A']*d0['W']
 
-    Inspecting the parameter estimates and the variance
+    Having setup our data, we can now define the psi function.
 
-    >>> mestimation.theta
-    >>> mestimation.variance
+    >>> def psi(theta):
+    >>>     return ee_aipw(theta,
+    >>>                    y=d['Y'],
+    >>>                    A=d['A'],
+    >>>                    W=d[['C', 'W']],
+    >>>                    X=d[['C', 'A', 'W', 'AW']],
+    >>>                    X1=d1[['C', 'A', 'W', 'AW']],
+    >>>                    X0=d0[['C', 'A', 'W', 'AW']])
 
-    More specifically, the average treatment effect and its variance are
+    Calling the M-estimation procedure. AIPW has 3 parameters with 2 coefficients in the propensity score model, and
+    4 coefficients in the outcome model, the total number of initial values should be 3+2+4=9. When Y is binary, it
+    will be best to start with ``[0., 0.5, 0.5, ...]`` followed by all ``0.`` for the initial values. Otherwise,
+    starting with all 0. as initials is reasonable.
 
-    >>> mestimation.theta[0]
-    >>> mestimation.variance[0, 0]
+    >>> estr = MEstimator(psi,
+    >>>                   init=[0., 0.5, 0.5, 0., 0., 0., 0., 0., 0.])
+    >>> estr.estimate(solver='lm')
 
-    The risk / mean had all been given treatment=1
+    Inspecting the parameter estimates, variance, and 95% confidence intervals
 
-    >>> mestimation.theta[1]
-    >>> mestimation.variance[1, 1]
+    >>> estr.theta
+    >>> estr.variance
+    >>> estr.confidence_intervals()
 
-    The risk / mean had all been given treatment=0
+    More specifically, the corresponding parameters are
 
-    >>> mestimation.theta[2]
-    >>> mestimation.variance[2, 2]
+    >>> estr.theta[0]     # causal mean difference of 1 versus 0
+    >>> estr.theta[1]     # causal mean under A=1
+    >>> estr.theta[2]     # causal mean under A=0
+    >>> estr.theta[3:5]   # propensity score regression coefficients
+    >>> estr.theta[5:]    # outcome regression coefficients
 
     References
     ----------
+    Hernán MA, & Robins JM. (2006). Estimating causal effects from epidemiological data.
+    *Journal of Epidemiology & Community Health*, 60(7), 578-586.
+
     Funk MJ, Westreich D, Wiesen C, Stürmer T, Brookhart MA, & Davidian M. (2011). Doubly robust estimation of causal
     effects. *American Journal of Epidemiology*, 173(7), 761-767.
 
     Tsiatis AA. (2006). Semiparametric theory and missing data. Springer, New York, NY.
     """
-    # Giving parameters short-hand names for simplicity of code
-    diff = theta[0]                 # This is the ATE
-    r1 = theta[1]                   # This is R1
-    r0 = theta[2]                   # This is R0
-    beta = theta[3:3+X.shape[1]]    # These are the regression coefficients for the outcome model
-    alpha = theta[3+X.shape[1]:]    # These are the regression coefficients for the treatment model
-
     # Ensuring correct typing
-    X = np.asarray(X)               # Convert to NumPy array
-    y = np.asarray(y)               # Convert to NumPy array
-
-    # Splitting X into A,W (treatment, covariates)
-    W = np.delete(X, treat_index, axis=1)             # Extract all-but treatment col A
-    A = X[:, treat_index].copy()                      # Extract treatment col A (MUST BE A COPY)
+    y = np.asarray(y)              # Convert to NumPy array
+    A = np.asarray(A)              # Convert to NumPy array
+    W = np.asarray(W)              # Convert to NumPy array
+    X = np.asarray(X)              # Convert to NumPy array
+    X1 = np.asarray(X1)            # Convert to NumPy array
+    X0 = np.asarray(X0)            # Convert to NumPy array
+
+    # Checking some shapes
+    if X.shape != X1.shape:
+        raise ValueError("The dimensions of X and X1 must be the same.")
+    if X.shape != X0.shape:
+        raise ValueError("The dimensions of X and X0 must be the same.")
+
+    # Extracting theta value for ease
+    mud = theta[0]                 # Parameter for average causal effect
+    mu1 = theta[1]                 # Parameter for the mean under A=1
+    mu0 = theta[2]                 # Parameter for the mean under A=0
+    alpha = theta[3:3+W.shape[1]]  # Parameter(s) for the propensity score model
+    beta = theta[3+W.shape[1]:]  # Parameter(s) for the outcome model
 
     # pi-model (logistic regression)
     pi_model = ee_logistic_regression(theta=alpha,    # Estimating logistic model
                                       X=W,
                                       y=A)
-    pi = inverse_logit(np.dot(W, alpha))              # Estimating the predicted probabilities for pi
+    pi = inverse_logit(np.dot(W, alpha))              # Estimating Pr(A|W)
+    if truncate is not None:                          # Truncating Pr(A|W) when requested
+        if truncate[0] > truncate[1]:
+            raise ValueError("truncate values must be specified in ascending order")
+        pi = np.clip(pi, a_min=truncate[0], a_max=truncate[1])
 
     # m-model (logistic regression)
-    if np.isin(y, [0, 1]).all() and not force_continuous:      # Checking outcome variable type is binary
-        m_model = ee_logistic_regression(theta=beta,           # Logistic regression parameters
-                                         X=X,
-                                         y=y)
-        X[:, treat_index] = 1                                  # Setting treatment to A=1 for all
-        ya1 = inverse_logit(np.dot(X, beta))                   # Calculating \hat{Y}(a=1)
-        X[:, treat_index] = 0                                  # Setting treatment to A=0 for all
-        ya0 = inverse_logit(np.dot(X, beta))                   # Calculating \hat{Y}(a=1)
+    # Checking outcome variable type
+    if np.isin(y, [0, 1]).all() and not force_continuous:
+        regression = ee_logistic_regression         # Use a logistic regression model
+        transform = inverse_logit                   # ... and need to inverse-logit transformation
     else:
-        m_model = ee_linear_regression(theta=beta,             # Linear regression parameters
-                                       X=X,
-                                       y=y)
-        X[:, treat_index] = 1                                  # Setting treatment to A=1 for all
-        ya1 = np.dot(X, beta)                                  # Calculating \hat{Y}(a=1)
-        X[:, treat_index] = 0                                  # Setting treatment A=0 for all
-        ya0 = np.dot(X, beta)                                  # Calculating \hat{Y}(a=0)
+        regression = ee_linear_regression           # Use a linear regression model
+        transform = identity                        # ... and need to apply the identity (no) transformation
+
+    m_model = regression(theta=beta,                # Estimating the outcome model
+                         y=y, X=X)
+    ya1 = transform(np.dot(X1, beta))               # Generating predicted values under X1
+    ya0 = transform(np.dot(X0, beta))               # Generating predicted values under X0
 
     # AIPW estimator
-    y1_star = (y*A/pi - ya1*(A-pi)/pi) - r1                    # Calculating \tilde{Y}(a=1)
-    y0_star = (y*(1-A)/(1-pi) + ya0*(A-pi)/(1-pi)) - r0        # Calculating \tilde{Y}(a=0)
-    ate = np.ones(y.shape[0]) * (r1 - r0) - diff               # Calculating the ATE
+    ace = np.ones(y.shape[0]) * (mu1 - mu0) - mud               # Calculating the ATE
+    y1_star = (y*A/pi - ya1*(A-pi)/pi) - mu1                    # Calculating \tilde{Y}(a=1)
+    y0_star = (y*(1-A)/(1-pi) + ya0*(A-pi)/(1-pi)) - mu0        # Calculating \tilde{Y}(a=0)
 
-    # Output (3+b+(b-1))-by-n array
-    return np.vstack((ate,               # theta[0] is for the ATE
+    # Output (3+b+c)-by-n array
+    return np.vstack((ace,               # theta[0] is for the ATE
                       y1_star[None, :],  # theta[1] is for R1
                       y0_star[None, :],  # theta[2] is for R0
-                      m_model,           # theta[3:b] is for the outcome model coefficients
-                      pi_model))         # theta[b:] is for the treatment model coefficients
+                      pi_model,          # theta[b] is for the treatment model coefficients
+                      m_model))          # theta[c] is for the outcome model coefficients
diff --git a/delicatessen/mestimation.py b/delicatessen/mestimation.py
index 23e0501..1155b56 100644
--- a/delicatessen/mestimation.py
+++ b/delicatessen/mestimation.py
@@ -1,14 +1,15 @@
 import numpy as np
 from scipy.optimize import newton, root
 from scipy.misc import derivative
+from scipy.stats import norm
 
 from delicatessen.utilities import partial_derivative
 
 
 class MEstimator:
-    r"""Generalized M-Estimator for stacked estimating equations.
+    r"""M-Estimator for stacked estimating equations.
 
-    M-estimation, or loosely referred to as estimating equations, is a general approach to point and variance
+    M-Estimation, or loosely referred to as estimating equations, is a general approach to point and variance
     estimation that consists of defining an estimator as the solution to an estimating equation (but does not require
     the derivative of a log-likelihood function). M-estimators satisify the following constraint
 
@@ -24,7 +25,7 @@ class MEstimator:
 
     M-Estimation consists of two broad step: point estimation and variance estimation. Point estimation is carried out
     by determining at which values of theta the given estimating equations are equal to zero. This is done via SciPy's
-    `newton` algorithm by default.
+    ``newton`` algorithm by default.
 
     For variance estimation, the sandwich asymptotic variance estimator is used, which consists of
 
@@ -64,8 +65,8 @@ class MEstimator:
     Parameters
     ----------
     stacked_equations : function, callable
-        Function that returns a b-by-n NumPy array of the estimating equations. See documentation for how to construct
-        a set of estimating equations.
+        Function that returns a b-by-n NumPy array of the estimating equations. See provided examples in the
+        documentation for how to construct a set of estimating equations.
     init : list, set, array
         Initial values to optimize for the function.
 
@@ -89,19 +90,20 @@ class MEstimator:
 
     Calling the M-estimation procedure
 
-    >>> mestimation = MEstimator(stacked_equations=psi, init=[0, 0, ])
-    >>> mestimation.estimate()
+    >>> estr = MEstimator(stacked_equations=psi, init=[0, 0, ])
+    >>> estr.estimate()
 
-    Inspecting the parameter estimates, the variance, and the asymptotic variance
+    Inspecting the parameter estimates, the variance, 95% confidence intervals, and the asymptotic variance
 
-    >>> mestimation.theta
-    >>> mestimation.variance
-    >>> mestimation.asymptotic_variance
+    >>> estr.theta
+    >>> estr.variance
+    >>> estr.confidence_intervals()
+    >>> estr.asymptotic_variance
 
     Alternatively, a custom estimating equation can be specified. This is done by constructing a valid estimating
-    equation for the ``MEstimator``. The ``MEstimator`` expects the `psi` function to return a b-by-n array, where b is
-    the number of parameters (length of theta) and n is the total number of observations. Below is an example of the
-    mean and variance estimating equation from before
+    equation for the ``MEstimator``. The ``MEstimator`` expects the ``psi`` function to return a b-by-n array, where b
+    is the number of parameters (length of ``theta``) and n is the total number of observations. Below is an example
+    of the mean and variance estimating equation from before
 
     >>> def psi(theta):
     >>>     y = np.array(y_dat)
@@ -111,12 +113,11 @@ class MEstimator:
 
     The M-estimation procedure is called using the same approach as before
 
-    >>> mestimation = MEstimator(stacked_equations=psi, init=[0, 0, ])
-    >>> mestimation.estimate()
+    >>> estr = MEstimator(stacked_equations=psi, init=[0, 0, ])
+    >>> estr.estimate()
 
     Note that ``len(init)`` should be equal to b. So in this case, two initial values are provided.
 
-    # TODO provide custom root-finding example here.
     Finally, the M-Estimator can also be run with a user-provided root-finding algorithm. To specify a custom
     root-finder, a function must be created by the user that consists of two keyword arguments (``stacked_equations``,
     ``init``) and must return only the optimized values. The following is an example with SciPy's Levenberg-Marquardt
@@ -132,8 +133,8 @@ class MEstimator:
     The provided custom root-finder can then be implemented like the following (continuing with the estimating equation
     from the previous example):
 
-    >>> mestimation = MEstimator(stacked_equations=psi, init=[0, 0, ])
-    >>> mestimation.estimate(solver=custom_solver)
+    >>> estr = MEstimator(stacked_equations=psi, init=[0, 0, ])
+    >>> estr.estimate(solver=custom_solver)
 
     References
     ----------
@@ -162,11 +163,12 @@ def estimate(self, solver='newton', maxiter=1000, tolerance=1e-9, dx=1e-9, order
         ----------
         solver : str, function, callable, optional
             Method to use for the root finding procedure. Default is the secant method (``scipy.optimize.newton``).
-            Another built-in option is the Levenberg-Marquardt algorithm (``scipy.optimize.root(method='lm')``).
-            Finally, any generic root-finding algorithm can be used via a user-provided callable object (function).
-            The function should consist of two keyword arguments: ``stacked_equations``, and ``init``. Additionally,
-            the function should return only the optimized values. Please review the example in the documentation or on
-            ReadTheDocs for how to provide a custom root-finding algorithm.
+            Other built-in option is the Levenberg-Marquardt algorithm (``scipy.optimize.root(method='lm')``), and
+            a modification of the Powell hybrid method (``scipy.optimize.root(method='hybr')``). Finally, any generic
+            root-finding algorithm can be used via a user-provided callable object (function). The function should
+            consist of two keyword arguments: ``stacked_equations``, and ``init``. Additionally, the function should
+            return only the optimized values. Please review the example in the documentation for how to provide a
+            custom root-finding algorithm.
         maxiter : int, optional
             Maximum iterations to consider for the root finding procedure. Default is 1000 iterations. For complex
             estimating equations (without preceding optimization), this value may need to be increased. This argument
@@ -177,15 +179,16 @@ def estimate(self, solver='newton', maxiter=1000, tolerance=1e-9, dx=1e-9, order
             this tolerance level (at this time). This argument is not used for user-specified solvers
         dx : float, optional
             Spacing to use to numerically approximate the partial derivatives of the bread matrix. Default is 1e-9,
-            which should work well for most applications. It is generally not recommended to have a large dx, since some
-            large values can poorly approximate derivatives.
+            which should work well for most applications. It is generally not recommended to have a large ``dx``, since
+            some large values can poorly approximate derivatives. Otherwise, also increase ``order``.
         order : int, optional
             Number of points to use to numerically approximate the partial derivative (must be an odd number). Default
             is 3, which is the default for SciPy.
         allow_pinv : bool, optional
             The default is ``True`` which uses ``numpy.linalg.pinv`` to find the inverse (or pseudo-inverse if matrix is
             non-invertible) for the bread. This default option is more robust to the possible matrices. If you want
-            to use ``numpy.linalg.inv`` instead (which does not support pseudo-inverse), set this parameter to False.
+            to use ``numpy.linalg.inv`` instead (which does not support pseudo-inverse), set this parameter to
+            ``False``.
 
         Returns
         -------
@@ -215,13 +218,10 @@ def estimate(self, solver='newton', maxiter=1000, tolerance=1e-9, dx=1e-9, order
         self.meat = np.dot(evald_theta, evald_theta.T) / self.n_obs         # Meat is a simple dot product of two arrays
 
         # Step 2.3: assembling the sandwich (variance)
-        if self.bread.ndim == 0:                        # NumPy's linalg throws an error if bread is a single value
-            bread_invert = 1 / self.bread               # ... so directly take inverse of the single value
-        else:                                           # otherwise bread must be a matrix
-            if allow_pinv:
-                bread_invert = np.linalg.pinv(self.bread)   # ... so find the inverse (or pseudo-inverse)
-            else:
-                bread_invert = np.linalg.inv(self.bread)    # ... so find the inverse (NOT pseudo-inverse)
+        if allow_pinv:                                  # re-worked to supports 1d theta
+            bread_invert = np.linalg.pinv(self.bread)   # ... so find the inverse (or pseudo-inverse)
+        else:
+            bread_invert = np.linalg.inv(self.bread)    # ... so find the inverse (NOT pseudo-inverse)
         # Two sets of matrix multiplication to get the sandwich variance
         sandwich = np.dot(np.dot(bread_invert, self.meat), bread_invert.T)
 
@@ -229,6 +229,43 @@ def estimate(self, solver='newton', maxiter=1000, tolerance=1e-9, dx=1e-9, order
         self.asymptotic_variance = sandwich       # Asymptotic variance estimate requires division by n (done above)
         self.variance = sandwich / self.n_obs     # Variance estimate requires division by n^2
 
+    def confidence_intervals(self, alpha=0.05):
+        r"""Calculate Wald-type :math:`(1 - \alpha) \times 100`% confidence intervals using the point estimates and
+        the sandwich variance. The formula for the confidence intervals are
+
+        .. math::
+
+            \hat{\theta} +/- Z_{\alpha / 2} \times \widehat{SE}(\hat{\theta})
+
+        Note
+        ----
+        The ``.estimate()`` function must be called before the confidence intervals can be calculated.
+
+        Parameters
+        ----------
+        alpha : float, optional
+            The :math:`\alpha` level for the corresponding confidence intervals. Default is 0.05, which calculate the
+            95% confidence intervals. Notice that :math:`0<\alpha<1`.
+
+        Returns
+        -------
+        array :
+            b-by-2 array, where row 1 is the confidence intervals for :math:`\theta_1`, ..., and row b is the confidence
+            intervals for :math:`\theta_b`
+        """
+        if not 0 < alpha < 1:
+            raise ValueError("`alpha` must be 0 < a < 1")
+
+        # 'Looking up' via Z table
+        z_alpha = norm.ppf(1 - alpha / 2, loc=0, scale=1)   # Z_alpha value for CI
+
+        # Calculating confidence intervals
+        param_se = np.sqrt(np.diag(self.variance))          # Take the diagonal of the sandwich and then SQRT
+        lower_ci = self.theta - z_alpha * param_se          # Calculate lower CI
+        upper_ci = self.theta + z_alpha * param_se          # Calculate upper CI
+
+        return np.asarray([lower_ci, upper_ci]).T
+
     def _mestimation_answer_(self, theta):
         """Internal function to evaluate the sum of the estimating equations. The summation must be internally evaluated
         since the bread requires calculation of the sum of the derivatives.
@@ -298,6 +335,16 @@ def _solve_coefficients_(stacked_equations, init, method, maxiter, tolerance):
             if opt.success == 0:
                 print("Root-finding failed to converge...")
                 raise RuntimeError(opt.message)
+        elif method in ['hybr', ]:
+            options = {"xtol": tolerance}
+            opt = root(stacked_equations,      # ... stacked equations to solve (should be written as sums)
+                       x0=np.asarray(init),    # ... initial values for solver
+                       method=method,          # ... allow for valid root-finding methods
+                       options=options)        # ... allow for options in hybrid
+            psi = opt.x                        # Error handling if fails to converge
+            if opt.success == 0:
+                print("Root-finding failed to converge...")
+                raise RuntimeError(opt.message)
         elif callable(method):
             try:
                 psi = method(stacked_equations=stacked_equations,
@@ -374,7 +421,7 @@ def _bread_matrix_(self, theta, stacked_equations, dx, order):
         # Evaluate the bread matrix
         if val_range == 1:                                       # When only a single theta is present
             d = derivative(stacked_equations, theta, dx=dx)      # ... approximate the derivative
-            return -1 * np.sum(d)                                # ... then return negative sum
+            return np.array([[-1 * np.sum(d), ], ])              # ... then return negative sum as 2d array
         else:                                                    # Otherwise approximate the partial derivatives
             bread_matrix = np.empty((val_range, val_range))      # ... create empty matrix
             for i in range(val_range):                           # ... for each i in len(theta)
diff --git a/delicatessen/utilities.py b/delicatessen/utilities.py
index 88bb84b..9393c91 100644
--- a/delicatessen/utilities.py
+++ b/delicatessen/utilities.py
@@ -26,7 +26,7 @@ def partial_derivative(func, var, point, output, dx, order):
     float
         Partial derivative at the variable at a particular points
     """
-    args = copy(point[:])    # Copy is needed here to prevent over-writing (not sure why SciPy overwrites if not)
+    args = copy(point[:])  # Copy is needed here to prevent over-writing (not sure why SciPy overwrites if not)
 
     def wraps(x):
         args[var] = x
@@ -36,11 +36,11 @@ def wraps(x):
 
 
 def logit(prob):
-    """Logit transformation of probabilities. Input can be a single probability of array of probabilities
+    """Logistic transformation of probabilities. Returns log-odds
 
     Parameters
     ----------
-    prob : float, array
+    prob : float, ndarray
         A single probability or an array of probabilities
 
     Returns
@@ -51,11 +51,11 @@ def logit(prob):
 
 
 def inverse_logit(logodds):
-    """Inverse logit transformation. Returns probabilities
+    """Inverse logistic transformation. Returns probabilities
 
     Parameters
     ----------
-    logodds : float, array
+    logodds : float, ndarray
         A single log-odd or an array of log-odds
 
     Returns
@@ -63,3 +63,23 @@ def inverse_logit(logodds):
     inverse-logit transformed results (i.e. probabilities for log-odds)
     """
     return 1 / (1 + np.exp(-logodds))
+
+
+def identity(value):
+    """Identity transformation. Returns itself
+
+    Note
+    ----
+    This function doesn't actually apply any transformation. It is used for arbitrary function calls that apply
+    transformations, and this is called when no transformation is to be applied
+
+    Parameters
+    ----------
+    value : float, ndarray
+        A single value or an array of values
+
+    Returns
+    -------
+    value
+    """
+    return value
diff --git a/delicatessen/version.py b/delicatessen/version.py
index 5ffcbf3..a4e2017 100644
--- a/delicatessen/version.py
+++ b/delicatessen/version.py
@@ -1,1 +1,1 @@
-__version__ = "0.1b3"
+__version__ = "0.1"
diff --git a/docs/Built-in Equations.rst b/docs/Built-in Equations.rst
index e1d6ca7..e8a279a 100644
--- a/docs/Built-in Equations.rst
+++ b/docs/Built-in Equations.rst
@@ -2,19 +2,20 @@ Built-in Estimating Equations
 '''''''''''''''''''''''''''''''''''''
 
 Here, we provide an overview of some of the built-in estimating equations with ``delicatessen``. This documentation is
-split into four sections, corresponding to basic, regression, survival, and causal estimating equations.
+split into three sections, corresponding to basic, regression, and causal estimating equations.
 
-All built-in estimating equations need to be 'wrapped' inside a outer function. Below is an example of an outer
-function, ``psi`` for the generic estimating equation ``ee`` (``ee`` does not exist but it a placeholder here).
+All built-in estimating equations need to be 'wrapped' inside an outer function. Below is a generic example of an outer
+function, where ``psi`` is the wrapper function and ``ee`` is the generic estimating equation example (``ee`` does not
+exist but it a placeholder here).
 
-.. code::python
+.. code::
 
     def psi(theta):
-        return ee(theta=theta, data=dat)
+        return ee(theta=theta, data=data)
 
 Here, the generic ``ee`` takes two inputs ``theta`` and ``data``. ``theta`` is the general theta vector that is present
 in all stacked estimating equations expected by ``delicatessen``. ``data`` is an argument that takes an input source
-of data. The ``dat`` provided should be **in the local scope** of the .py file this function lives in.
+of data. The ``data`` provided should be **in the local scope** of the ``.py`` file this function lives in.
 
 After wrapped in an outer function, the function can be passed to ``MEstimator``. See the examples below for further
 details and examples.
@@ -22,19 +23,20 @@ details and examples.
 Basic Equations
 =============================
 
-Some basic estimating equations for the mean and variance are provided.
+Some basic estimating equations are provided.
 
 Mean
 ----------------------------
 
-The most basic is the estimating equation for the mean: ``ee_mean``. To illustrate, consider we wanted to estimated the
-mean for the following data
+The most basic available estimating equation is for the mean: ``ee_mean``. To illustrate, consider we wanted to
+estimated the mean for the following data
 
 .. code::
 
     obs_vals = [1, 2, 1, 4, 1, 4, 2, 4, 2, 3]
 
-To use ``ee_mean`` with ``MEstimator``, this function will be wrapped in an outer function. Below is an example
+To use ``ee_mean`` with ``MEstimator``, this function will be wrapped in an outer function. Below is an illustration of
+this wrapper function
 
 .. code::
 
@@ -51,17 +53,20 @@ After creating the wrapper function, the M-Estimator can be called like the foll
 
     from delicatessen import MEstimator
 
-    mestimation = MEstimator(stacked_equations=psi, init=[0, ])
-    mestimation.estimate()
+    estr = MEstimator(stacked_equations=psi, init=[0, ])
+    estr.estimate()
+
+    print(estr.theta)   # [2.4, ])
 
 Since ``ee_mean`` consists of a single parameter, only a single ``init`` value is provided.
 
 Robust Mean
 ----------------------------
 
-Sometimes extreme outliers are observed. The mean will be sensitive to these extreme outliers, but excluding them also
-seems like 'cheating'. Instead, a robust version of the mean could be considered. Consider the following generic data,
-where there are two extreme outliers
+Sometimes extreme observations, termed outliers, occur. The mean is generally sensitive to these outliers. A common
+approach to handling outliers is to exclude them. However, exclusion ignores all information contributed by outliers,
+and should only be done when outliers are a result of experimental error. Robust statistics have been proposed as
+middle ground, whereby outliers contribute to estimation but their influence is constrained.
 
 .. code::
 
@@ -71,8 +76,7 @@ Rather than excluding the -10 and 12, we can use the robust mean proposed by Hub
 pre-specified level. Therefore, they still contribute information but only values up to the bound. In this example, a
 bound of -6,6 will be applied.
 
-The robust mean estimating equation is available in ``ee_mean_robust``. Below is an example (including the wrapper
-function and call).
+The robust mean estimating equation is available in ``ee_mean_robust``.
 
 .. code::
 
@@ -82,13 +86,11 @@ function and call).
     def psi(theta):
         return ee_mean_robust(theta=theta, y=obs_vals, k=6)
 
-    mestimation = MEstimator(stacked_equations=psi, init=[0, ])
-    mestimation.estimate()
+    estr = MEstimator(stacked_equations=psi, init=[0, ])
+    estr.estimate()
 
-    print("Mean:    ", mestimation.theta)
-    print("Variance:", mestimation.variance)
+    print(estr.theta)  # [2.0, ]
 
-Therefore, the robust mean point and variance estimate are displayed.
 
 Mean and Variance
 ----------------------------
@@ -110,11 +112,10 @@ The mean-variance estimating equation can be implemented as follows (remember th
     def psi(theta):
         return ee_mean_variance(theta=theta, y=obs_vals)
 
-    mestimation = MEstimator(stacked_equations=psi, init=[0, 1, ])
-    mestimation.estimate()
+    estr = MEstimator(stacked_equations=psi, init=[0, 1, ])
+    estr.estimate()
 
-    print("theta:     ", mestimation.theta)
-    print("Var(theta):", mestimation.variance)
+    print(estr.theta)  # [2.4, 1.44]
 
 *Note* the ``init`` here takes two values because the stacked estimating equations has a length of 2 (``theta`` is
 b-by-1 where b=2). The first value of ``theta`` is the mean and the second is the variance. Now, the variance output
@@ -124,14 +125,12 @@ estimated variance of the mean and the second is the estimated variance of the v
 Regression
 =============================
 
-Several basic regression model estimating equations are provided.
+Several common regression models are provided as built-in estimating equations.
 
 Linear Regression
 ----------------------------
 
-The estimating equations for linear regression predict a continuous outcome as a function of provided covariates. The
-implementation of linear regression here is similar to ordinary least squares, but the variance here is robust.
-Specifically, the sandwich variance estimator of M-Estimation is robust.
+The estimating equations for linear regression predict a continuous outcome as a function of provided covariates.
 
 To demonstrate application, consider the following simulated data set
 
@@ -147,9 +146,9 @@ To demonstrate application, consider the following simulated data set
     data['Y'] = 0.5 + 2*data['X'] - 1*data['Z'] + np.random.normal(loc=0, size=n)
     data['C'] = 1
 
-In this case, X and Z are the independent variables and Y is the dependent variable. Here C is necessary as a column
-since we need to manually provide the intercept (this may be different from other formula-based packages that
-automatically add the intercept to the regression).
+In this case, ``X`` and ``Z`` are the independent variables and ``Y`` is the dependent variable. Here the column ``C``
+is created to be the intercept column, since the intercept needs to be manually provided (this may be different from
+other formula-based packages that automatically add the intercept to the regression).
 
 For this data, we can now create the wrapper function for the ``ee_linear_regression`` estimating equations
 
@@ -168,15 +167,15 @@ and their variance
 
 .. code::
 
-    mestimation = MEstimator(stacked_equations=psi, init=[0., 0., 0.])
-    mestimation.estimate()
+    estr = MEstimator(stacked_equations=psi, init=[0., 0., 0.])
+    estr.estimate()
 
-    print("theta:     ", mestimation.theta)
-    print("Var(theta):", mestimation.variance)
+    print(estr.theta)
+    print(estr.variance)
 
-Note that ``X`` is 3 covariates, meaning ``init`` needs 3 starting values. The linear regression done here should match
-the ``statsmodels`` generalized linear model with a robust variance estimate. Below is code demonstrating how to
-estimate the same quantities with ``statsmodels.glm``.
+Note that there are 3 independent variables, meaning ``init`` needs 3 starting values. The linear regression done here
+should match the ``statsmodels`` generalized linear model with a robust variance estimate. Below is code on how to
+compare to ``statsmodels.glm``.
 
 .. code::
 
@@ -194,10 +193,8 @@ causal section.
 Logistic Regression
 ----------------------------
 
-In the case of a binary dependent variable, logistic regression can instead be performed (no linear probability models
-here!).
-
-To demonstrate application, consider the following simulated data set
+In the case of a binary dependent variable, logistic regression can instead be performed. Consider the following
+simulated data set
 
 .. code::
 
@@ -212,9 +209,9 @@ To demonstrate application, consider the following simulated data set
     data['Y'] = np.random.binomial(n=1, p=logistic.cdf(0.5 + 2*data['X'] - 1*data['Z']), size=n)
     data['C'] = 1
 
-In this case, X and Z are the independent variables and Y is the dependent variable. Here C is necessary as a column
-since we need to manually provide the intercept (this may be different from other formula-based packages that
-automatically add the intercept to the regression).
+In this case, ``X`` and ``Z`` are the independent variables and ``Y`` is the dependent variable. Here the column ``C``
+is created to be the intercept column, since the intercept needs to be manually provided (this may be different from
+other formula-based packages that automatically add the intercept to the regression).
 
 For this data, we can now create the wrapper function for the ``ee_logistic_regression`` estimating equations
 
@@ -233,15 +230,15 @@ and their variance
 
 .. code::
 
-    mestimation = MEstimator(stacked_equations=psi, init=[0., 0., 0.])
-    mestimation.estimate()
+    estr = MEstimator(stacked_equations=psi, init=[0., 0., 0.])
+    estr.estimate()
 
-    print("theta:     ", mestimation.theta)
-    print("Var(theta):", mestimation.variance)
+    print(estr.theta)
+    print(estr.variance)
 
-Note that ``X`` is 3 covariates, meaning ``init`` needs 3 starting values. The logistic regression done here should
-match the ``statsmodels`` generalized linear model with a robust variance estimate. Below is code demonstrating how to
-estimate the same quantities with ``statsmodels.glm``.
+Note that there are 3 independent variables, meaning ``init`` needs 3 starting values. The logistic regression done here
+should match the ``statsmodels`` generalized linear model with a robust variance estimate. Below is code on how to
+compare to ``statsmodels.glm``.
 
 .. code::
 
@@ -258,16 +255,236 @@ equations can be stacked together (including multiple regression models). This a
 causal section.
 
 
+Dose-Response
+=============================
+
+Estimating equations for dose-response relationships are also included. The following examples use the data from
+Inderjit et al. (2002). This data can be loaded via
+
+.. code::
+
+    d = load_inderjit()   # Loading array of data
+    dose_data = d[:, 1]   # Dose data
+    resp_data = d[:, 0]   # Response data
+
+
+4-parameter Logistic
+----------------------------
+
+The 4-parameter logistic model (4PL) consists of parameters for the lower-limit of the response, the effective dose,
+steepness of the curve, and the upper-limit of the response.
+
+The wrapper function for the 4PL model should look like
+
+.. code::
+
+    from delicatessen import MEstimator
+    from delicatessen.estimating_equations import ee_4p_logistic
+
+    def psi(theta):
+        # Estimating equations for the 4PL model
+        return ee_4p_logistic(theta=theta, X=dose_data, y=resp_data)
+
+After creating the wrapper function, we can now call the M-Estimation procedure to estimate the coefficients for the
+4PL model and their variance
+
+.. code::
+
+    estr = MEstimator(psi, init=[np.min(resp_data),
+                                 (np.max(resp_data)+np.min(resp_data)) / 2,
+                                 (np.max(resp_data)+np.min(resp_data)) / 2,
+                                 np.max(resp_data)])
+    estr.estimate(solver='lm')
+
+    print(estr.theta)
+    print(estr.variance)
+
+When you use 4PL, you may notice convergence errors. This estimating equation can be hard to optimize since it has
+implicit bounds the root-finder isn't aware of. To avoid these issues, we can give the root-finder good starting values.
+
+First, the upper limit should *always* be greater than the lower limit. Second, the ED50 should be between the lower
+and upper limits. Third, the sign for the steepness depends on whether the response declines (positive) or the response
+increases (negative). Finally, some solvers may be better suited to the problem, so try a few different options. With
+decent initial values, I have found ``lm`` to be fairly reliable.
+
+For the 4PL, good general starting values I have found are the following. For the lower-bound, give the minimum response
+value as the initial. For ED50, give the mid-point between the maximum response and the minimum response. The initial
+value for steepness is more difficult. Ideally, we would give a starting value of zero, but that will fail in this
+4PL. Giving a larger starting value (between 2 to 8) works in this example. For the upper-bound, give the maximum
+response value as the initial.
+
+To summarize, be sure to examine your data (e.g., scatterplot). This will help to determine the initial starting values
+for the root-finding procedure. Otherwise, you may come across a convergence error.
+
+
+3-parameter Logistic
+----------------------------
+
+The 3-parameter logistic model (3PL) consists of parameters for the effective dose, steepness of the curve, and the
+upper-limit of the response. Here, the lower-limit is pre-specified and is no longer being estimated.
+
+The wrapper function for the 3PL model should look like
+
+.. code::
+
+    from delicatessen import MEstimator
+    from delicatessen.estimating_equations import ee_3p_logistic
+
+    def psi(theta):
+        # Estimating equations for the 3PL model
+        return ee_3p_logistic(theta=theta, X=dose_data, y=resp_data,
+                              lower=0)
+
+Since the shortest a root of a plant could be zero, a lower limit of zero makes sense here.
+
+After creating the wrapper function, we can now call the M-Estimation procedure to estimate the coefficients for the
+3PL model and their variance
+
+.. code::
+
+    estr = MEstimator(psi, init=[(np.max(resp_data)+np.min(resp_data)) / 2,
+                                 (np.max(resp_data)+np.min(resp_data)) / 2,
+                                 np.max(resp_data)])
+    estr.estimate(solver='lm')
+
+    print(estr.theta)
+    print(estr.variance)
+
+As before, you may notice convergence errors. This estimating equation can be hard to optimize since it has implicit
+bounds the root-finder isn't aware of. To avoid these issues, we can give the root-finder good starting values.
+
+For the 3PL, good general starting values I have found are the following. For ED50, give the mid-point between the
+maximum response and the minimum response. The initial value for steepness is more difficult. Ideally, we would give a
+starting value of zero, but that will fail in this 3PL. Giving a larger starting value (between 2 to 8) works in this
+example. For the upper-bound, give the maximum response value as the initial.
+
+To summarize, be sure to examine your data (e.g., scatterplot). This will help to determine the initial starting values
+for the root-finding procedure. Otherwise, you may come across a convergence error.
+
+2-parameter Logistic
+----------------------------
+
+The 2-parameter logistic model (2PL) consists of parameters for the effective dose, and steepness of the curve. Here,
+the lower-limit and upper-limit are pre-specified and no longer being estimated.
+
+The wrapper function for the 3PL model should look like
+
+.. code::
+
+    from delicatessen import MEstimator
+    from delicatessen.estimating_equations import ee_2p_logistic
+
+    def psi(theta):
+        # Estimating equations for the 2PL model
+        return ee_2p_logistic(theta=theta, X=dose_data, y=resp_data,
+                              lower=0, upper=8)
+
+While a lower-limit of zero makes sense in this example, the upper-limit of 8 is poorly motivated (and thus this should
+only be viewed as an example of the 2PL model and not how it should be applied in practice). Setting the limits as
+constants should be motivated by substantive knowledge of the problem.
+
+After creating the wrapper function, we can now call the M-Estimation procedure to estimate the coefficients for the
+2PL model and their variance
+
+.. code::
+
+    estr = MEstimator(psi, init=[(np.max(resp_data)+np.min(resp_data)) / 2,
+                                 (np.max(resp_data)+np.min(resp_data)) / 2])
+    estr.estimate(solver='lm')
+
+    print(estr.theta)
+    print(estr.variance)
+
+As before, you may notice convergence errors. This estimating equation can be hard to optimize since it has implicit
+bounds the root-finder isn't aware of. To avoid these issues, we can give the root-finder good starting values.
+
+For the 2PL, good general starting values I have found are the following. For ED50, give the mid-point between the
+maximum response and the minimum response. The initial value for steepness is more difficult. Ideally, we would give a
+starting value of zero, but that will fail in this 2PL.
+
+To summarize, be sure to examine your data (e.g., scatterplot). This will help to determine the initial starting values
+for the root-finding procedure. Otherwise, you may come across a convergence error.
+
+
+ED(:math:`\delta`)
+----------------------------
+
+In addition to the :math:`x`-parameter logistic models, an estimating equation to estimate a corresponding
+:math:`\delta` effective dose is available. Notice that this estimating equation should be stacked with one of
+the :math:`x`-PL models. Here, we demonstrate with the 3PL model.
+
+Here, our interest is in the following effective doses: 0.05, 0.10, 0.20, 0.80. The wrapper function for the 3PL model
+and estimating equations for these effective doses are
+
+.. code::
+
+    def psi(theta):
+        lower_limit = 0
+
+        # Estimating equations for the 3PL model
+        pl3 = ee_3p_logistic(theta=theta, X=d[:, 1], y=d[:, 0],
+                             lower=lower_limit)
+
+        # Estimating equations for the effective concentrations
+        ed05 = ee_effective_dose_delta(theta[3], y=resp_data, delta=0.05,
+                                       steepness=theta[0], ed50=theta[1],
+                                       lower=lower_limit, upper=theta[2])
+        ed10 = ee_effective_dose_delta(theta[4], y=resp_data, delta=0.10,
+                                       steepness=theta[0], ed50=theta[1],
+                                       lower=lower_limit, upper=theta[2])
+        ed20 = ee_effective_dose_delta(theta[5], y=resp_data, delta=0.20,
+                                       steepness=theta[0], ed50=theta[1],
+                                       lower=lower_limit, upper=theta[2])
+        ed80 = ee_effective_dose_delta(theta[6], y=resp_data, delta=0.80,
+                                       steepness=theta[0], ed50=theta[1],
+                                       lower=lower_limit, upper=theta[2])
+
+        # Returning stacked estimating equations
+        return np.vstack((pl3,
+                          ed05,
+                          ed10,
+                          ed20,
+                          ed80))
+
+Notice that the estimating equations are stacked together in the order of the ``theta`` vector.
+
+After creating the wrapper function, we can now call the M-Estimation procedure to estimate the coefficients for the
+3PL model, the ED for the :math:`\delta` values, and their variance
+
+.. code::
+
+    midpoint = (np.max(resp_data)+np.min(resp_data)) / 2
+    estr = MEstimator(psi, init=[midpoint,
+                                 midpoint,
+                                 np.max(resp_data),
+                                 midpoint,
+                                 midpoint,
+                                 midpoint,
+                                 midpoint])
+    estr.estimate(solver='lm')
+    print(estr.theta)
+    print(estr.variance)
+
+Since the ED for :math:`\delta`'s are transformations of the other parameters, there starting values are less important
+(the root-finders are better at solving those equations). Again, we can make it easy on the solver by having the
+starting point for each being the mid-point of the response values.
+
+
 Causal Inference
 =============================
 
-To demonstrate the utility of M-estimation, particularly how estimating equations can be 'stacked' together, then
-still have an appropriate variance estimator, several causal inference estimators are provided here.
+This next section describes a the available estimators for the causal mean. These estimators all rely on specific
+identification conditions to be able to interpret the estimate of the mean (or mean difference) as an estimate of the
+causal mean. For information on these assumptions, I recommend this
+`paper<https://www.ncbi.nlm.nih.gov/labs/pmc/articles/PMC2652882/>`_ as a general introduction.
 
-It is recommended that you are familiar with causal inference (particularly the identification conditions of these
-estimators) before using this utility widely. Causal inference is a difficult endeavour, my dear user!
+This section procedures that the identification conditions have been previously deliberated, and the causal mean is
+identified and is estimable (see this `paper<https://arxiv.org/abs/2108.11342>`_ or this
+`paper<https://arxiv.org/abs/1904.02826>`_ for more information on this concept).
 
-In the following examples, we will use the generic data example here, where W is a confounder of the A-Y relationship
+With that aside, let's proceed through the available estimators of the causal means. In the following examples, we will
+use the generic data example here, where :math:`Y(a)` is independent of :math:`A` conditional on :math:`W`. Below is
+a sample data set
 
 .. code::
 
@@ -280,231 +497,372 @@ In the following examples, we will use the generic data example here, where W is
     d['Y'] = (1-d['A'])*d['Ya0'] + d['A']*d['Ya1']
     d['C'] = 1
 
-Now to the examples
+Here, we don't get to see the potential outcomes :math:`Y(a)`, but instead estimate the mean under different plans
+using the observed data, :math:`Y,A,W`.
 
-G-computation
-----------------------------
+Inverse probability weighting
+-------------------------------------
 
-First, is g-computation. The built-in estimating equations for g-computation calculate the average treatment effect,
-risk / mean under all-treated, and the risk / mean under none-treated.
+First, we use the inverse probability weighting (IPW) estimator, which models the probability of :math:`A` conditional
+on :math:`W`. In general, the IPW estimator for the mean difference can be written as
 
-*A limitation*: the g-computation, as implemented in the built-in estimating equation only uses a single outcome model
-and that outcome model does *not* support interaction terms. Here the g-computation is meant as a basic example. For
-more general use, the provided estimating equation should be adapted. But the built-in estimating equation will provide
-a basic structure for user's to build off of.
+.. math::
 
-To load the estimating equations, we call
+    \frac{1}{n} \sum_{i=1}^n \frac{Y_i A_i}{Pr(A_i = 1 | W_i; \hat{\alpha})} - \frac{1}{n}
+    \sum_{i=1}^n \frac{Y_i (1-A_i)}{Pr(A_i = 0 | W_i; \hat{\alpha})}
+
+In ``delicatessen``, the built-in IPW estimator consists of 4 estimating equations, with both binary and continuous
+outcomes supported by ``ee_ipw`` (since we are using the Horwitz-Thompson estimator). The stacked estimating equations
+are
+
+.. math::
+
+    \sum_i^n \psi_d(Y_i, A_i, \pi_i, \theta_0) = \sum_i^n (\theta_1 - \theta_2) - \theta_0 = 0
+
+    \sum_i^n \psi_1(Y_i, A_i, \pi_i, \theta_1) = \sum_i^n \frac{A_i \times Y_i}{\pi_i} - \theta_1 = 0
+
+    \sum_i^n \psi_0(Y_i, A_i, \pi_i, \theta_2) = \sum_i^n \frac{(1-A_i) \times Y_i}{1-\pi_i} - \theta_2 = 0
+
+    \sum_i^n \psi_g(A_i, W_i, \theta) = \sum_i^n (A_i - expit(W_i^T \alpha)) W_i = 0
+
+where :math:`\theta_1` is the average causal effect, :math:`\theta_2` is the mean under the plan where
+:math:`A=1` for everyone, :math:`\theta_3` is the mean under the plan where :math:`A=0` for everyone, and
+:math:`\alpha` is the parameters for the logistic model used to estimate the propensity scores.
+
+To load the estimating equations,
 
 .. code::
 
     from delicatessen import MEstimators
-    from delicatessen.estimating_equations import ee_gformula
+    from delicatessen.estimating_equations import ee_ipw
 
-Again, we will wrap the built-in estimating equations inside a function.
+The estimating equation is then wrapped inside the wrapper ``psi`` function. Notice that the estimating equation has
+4 non-optional inputs: the parameter values, the outcomes, the actions, and the covariates to model the propensity
+scores with.
 
 .. code::
 
     def psi(theta):
-        return ee_gformula(theta, X=d[['C', 'A', 'W']], y=d['Y'], treat_index=1)
+        return ee_ipw(theta,                 # Parameters
+                      y=d['Y'],              # Outcome
+                      A=d['A'],              # Action (exposure, treatment, etc.)
+                      W=d[['C', 'W']])       # Covariates for PS model
 
-The arguments for ``ee_gformula`` are the :math:`\theta` values, the covariates (including an intercept (C) and the
-treatment (A)), the outcome values (Y), and the column index for the treatment in X. Here, 1 designates the second
-column (python uses zero-indexing), which corresponds to 'A' in how the X data is formatted.
+Note that we add an intercept to the logistic model by adding a column of 1's via ``d['C']``.
 
-Now we can call the M-estimator to solve for the values and the variance. Here, the initial values provided must be
-3+*b* (where *b* is the number of columns in X). This is because the g-computation estimating equations output the
-average treatment effect, risk under all-treated, risk under none-treated, and the regression model coefficients.
+Here, the initial values provided must be 3+*b* (where *b* is the number of columns in W). For binary
+outcomes, it will likely be best practice to have the initial values set as ``[0., 0.5, 0.5, ...]``. followed by b
+``0.``'s. For continuous outcomes, all ``0.`` can be used instead. Furthermore, a logistic model for the propensity
+scores could be optimized outside of ``delicatessen`` and those (pre-washed) regression estimates can be passed as
+initial values to speed up optimization.
 
-As for starting values, it will likely be best practice to have the initial values set as  [0., 0.5, 0.5, ...] in
-general. The regression initial values can also be pre-washed to speed up optimization.
+Now we can call the M-estimator to solve for the values and the variance.
 
 .. code::
 
-    mestimation = MEstimator(stacked_equations=psi, init=[0., 0.5, 0.5, 0., 0., 0.])
-    mestimation.estimate(solver='lm')
+    estr = MEstimator(psi, init=[0., 0.5, 0.5, 0., 0.])
+    estr.estimate(solver='lm')
 
-Now the average treatment effect, as well as the variance, can be output. Here, a key advantage of M-estimation can be
-seen. The form of an M-estimator allows us to estimate the variance directly, while appropriately allowing for the
-uncertainty in the regression model parameters to be carried forward. M-estimation does this automatically for us.
-Essentially, we do not need to bootstrap to estimate the variance!
+After successful optimization, we can inspect the estimated values.
 
 .. code::
 
-    mestimation.theta[0]
-    mestimation.variance[0, 0]
+    estr.theta[0]    # causal mean difference of 1 versus 0
+    estr.theta[1]    # causal mean under X1
+    estr.theta[2]    # causal mean under X0
+    estr.theta[3:]   # logistic regression coefficients
 
-Besides the average treatment effect, the risk / mean under all-treated can be extracted by
+The variance and Wald-type confidence intervals can also be output via
 
 .. code::
 
-    mestimation.theta[1]
-    mestimation.variance[1, 1]
+    estr.variance
+    estr.confidence_intervals()
 
-and the risk / mean under none-treated can be extracted by
+The IPW estimators demonstrates a key advantage of M-Estimation. The stacked estimating equations means that the
+sandwich variance correctly incorporates the uncertainty in estimation of the propensity scores into the parameter(s)
+of interest (e.g., average causal effect). Therefore, we do not have to rely on the nonparametric bootstrap
+(computationally cumbersome) or the GEE-trick (conservative estimate of the variance for the average causal effect).
 
-.. code::
 
-    mestimation.theta[2]
-    mestimation.variance[2, 2]
+G-computation
+----------------------------
 
-The ``ee_gformula`` supports both binary and continuous outcomes. Inside the function, it automatically detects whether
-the outcome data is binary. If the outcome data is not binary, then it defaults to using a linear regression model
-(but you can also force the use of a linear regression model for binary data by setting ``force_continuous=True``
+Second, we use g-computation, which models :math:`Y` conditional on :math:`A` and :math:`W`. In general, g-computation
+for the mean difference can be written as
 
-To summarize, the key advantage of M-estimation here is that it *appropriately* estimates the variance. We do *not*
-need to bootstrap in this case (and more generally if the sample size is sufficiently large).
+.. math::
 
-Inverse probability weighting
--------------------------------------
+    \frac{1}{n} \sum_{i=1}^n m_1(W_i; \hat{\beta}) - \frac{1}{n} \sum_{i=1}^n m_0(W_i; \hat{\beta})
+
+where :math:`m_a(W_i; \hat{\beta}) = E[Y_i|A_i=a,W_i; \hat{\beta}]`. In ``delicatessen``, the built-in g-computation
+consists of either 2 estimating equations or 4 estimating equations, with both binary and continuous outcomes supported.
+The 2 stacked estimating equations are
+
+.. math::
+
+    \sum_i^n \psi_1(Y_i, X_i, \theta_1) = \sum_i^n g(\hat{Y}_i^a) - \theta_1 = 0
+
+    \sum_i^n \psi_m(Y_i, X_i, \theta) = \sum_i^n (Y_i - \text{expit}(X_i^T \theta)) X_i = 0
+
+
+where :math:`\theta_1` is the mean under the action :math:`a`, and :math:`\beta` is the parameters for the regression
+model used to estimate the outcomes. Notice that the g-computation procedure supports generic deterministic plans
+(e.g., set :math:`A=1` for all, set :math:`A=0` for all, set :math:`A=1` if :math:`W=1` otherwise :math:`A=0`, etc.).
+These plans are more general than those allowed by either the built-in IPW or built-in AIPW estimating equations.
+
+The 4 stacked estimating equations instead compare the mean difference between two action plans. The estimating
+equations are
+
+.. math::
+
+    \sum_i^n \psi_1(Y_i, X_i, \theta_1) = \sum_i^n (\theta_2 - \theta_3) - \theta_1 = 0
+
+    \sum_i^n \psi_1(Y_i, X_i, \theta_2) = \sum_i^n g(\hat{Y}_i^a) - \theta_2 = 0
+
+    \sum_i^n \psi_1(Y_i, X_i, \theta_3) = \sum_i^n g(\hat{Y}_i^a) - \theta_3 = 0
+
+    \sum_i^n \psi_m(Y_i, X_i, \theta) = \sum_i^n (Y_i - \text{expit}(X_i^T \theta)) X_i = 0
 
-Rather than modeling the outcome, we can choose the inverse probability weighting (IPW) estimator, which models the
-probability of treatment. The estimating equations for the IPW estimator are also built-in to ``delicatessen``.
 
-To load the estimating equations, we call
+where :math:`\theta_1` is the average causal effect, :math:`\theta_2` is the mean under the first plan, :math:`\theta_3`
+is the mean under the second, and :math:`\beta` is the parameters for the regression model used to predict the
+outcomes.
+
+To load the estimating equations,
 
 .. code::
 
     from delicatessen import MEstimators
-    from delicatessen.estimating_equations import ee_ipw
+    from delicatessen.estimating_equations import ee_gformula
 
-As with every built-in estimating equation, we will wrap it inside a function.
+The estimating equation is then wrapped inside the wrapper ``psi`` function. In the first example, we focus on
+estimating the average causal effect. Notice that for ``ee_gformula`` some additional data prep is necessary.
+Specifically, we need to create a copy of the data set where ``A`` is set to the value our plan dictates
+(e.g., ``A=1``). Below is code that does this step and creates the wrapper function
 
 .. code::
 
+    # Creating data under the plans
+    d1 = d.copy()
+    d1['A'] = 1
+    d0 = d.copy()
+    d0['A'] = 0
+
+    # Creating interaction terms
+    d['AxW'] = d['A'] * d['W']
+    d1['AxW'] = d1['A'] * d1['W']
+    d0['AxW'] = d0['A'] * d0['W']
+
     def psi(theta):
-        return ee_ipw(theta, X=d[['C', 'A', 'W']], y=d['Y'], treat_index=1)
+        return ee_gformula(theta,                        # Parameters
+                           y=d['Y'],                     # Outcome
+                           X=d[['C', 'A', 'W', 'AxW']],  # Observed
+                           X=d1[['C', 'A', 'W', 'AxW']], # Plan 1
+                           X=d0[['C', 'A', 'W', 'AxW']]) # Plan 2
 
-The arguments for ``ee_ipw`` are the :math:`\theta` values, the covariates (including an intercept (C) and the treatment
-(A)), the outcome values (Y), and the column index for the treatment in X. Here, 1 designates the second column
-(python uses zero-indexing), which corresponds to 'A' in how the X data is formatted.
+Note that we add an intercept to the outcome model by adding a column of 1's via ``d['C']``.
+
+Here, the initial values provided must be 3+*b* (where *b* is the number of columns in X). For binary
+outcomes, it will likely be best practice to have the initial values set as ``[0., 0.5, 0.5, ...]``. followed by b
+``0.``'s. For continuous outcomes, all ``0.`` can be used instead. Furthermore, a regression model for the outcomes
+could be optimized outside of ``delicatessen`` and those (pre-washed) regression estimates can be passed as
+initial values to speed up optimization.
+
+Now we can call the M-estimator to solve for the values and the variance.
+
+.. code::
 
-Now we can call the M-estimator to solve for the values and the variance. Here, the initial values provided must be
-3+*b* (where *b* is the number of columns in X *minus 1*). This is because the IPW estimating equations output the
-average treatment effect, risk under all-treated, risk under none-treated, and the logistic regression model
-coefficients. Since we are modeling the conditional probability of A, one column in X is excluded from the covariates.
+    estr = MEstimator(psi, init=[0., 0.5, 0.5, 0., 0., 0., 0.])
+    estr.estimate(solver='lm')
 
-As for starting values, it will likely be best practice to have the initial values set as  [0., 0.5, 0.5, ...] in
-general. The regression initial values can also be pre-washed to speed up optimization.
+After successful optimization, we can inspect the estimated values.
 
 .. code::
 
-    mestimation = MEstimator(stacked_equations=psi, init=[0., 0.5, 0.5, 0., 0., 0.])
-    mestimation.estimate(solver='lm')
+    estr.theta[0]    # causal mean difference of 1 versus 0
+    estr.theta[1]    # causal mean under X1
+    estr.theta[2]    # causal mean under X0
+    estr.theta[3:]   # regression coefficients
 
-Now the average treatment effect, as well as the variance, can be output. Here, a key advantage of M-estimation can be
-seen. The form of an M-estimator allows us to estimate the variance directly, while appropriately allowing for the
-uncertainty in the regression model parameters to be carried forward. M-estimation does this automatically for us.
-Essentially, we do not need to bootstrap or use the GEE-trick for IPW to estimate the variance!
+The variance and Wald-type confidence intervals can also be output via
 
 .. code::
 
-    mestimation.theta[0]
-    mestimation.variance[0, 0]
+    estr.variance
+    estr.confidence_intervals()
 
-Besides the average treatment effect, the risk / mean under all-treated can be extracted by
+Again, a key advantage of M-Estimation is demonstrated here. The stacked estimating equations means that the
+sandwich variance correctly incorporates the uncertainty in estimation of the outcome model into the parameter(s)
+of interest (e.g., average causal effect). Therefore, we do not have to rely on the nonparametric bootstrap
+(computationally cumbersome).
+
+As a second example, we now demonstrate the flexbility of ``ee_gformula`` to estimate other plans. Here, we estimate
+the causal mean under the plan where only those with :math:`W=1` have :math:`A=1`. As before, we need to generate
+this distribution of covariates and wrap the built-in estimating equations.
 
 .. code::
 
-    mestimation.theta[1]
-    mestimation.variance[1, 1]
+    # Creating data under the plans
+    da = d.copy()
+    da['A'] = np.where(da['W'] == 1, 1, 0)
+
+    # Creating interaction terms
+    d['AxW'] = d['A'] * d['W']
+    da['AxW'] = da['A'] * da['W']
 
-and the risk / mean under none-treated can be extracted by
+    def psi(theta):
+        return ee_gformula(theta,                        # Parameters
+                           y=d['Y'],                     # Outcome
+                           X=d[['C', 'A', 'W', 'AxW']],  # Observed
+                           X=da[['C', 'A', 'W', 'AxW']]) # Plan
+
+Now we can call the M-estimator to solve for the values and the variance.
 
 .. code::
 
-    mestimation.theta[2]
-    mestimation.variance[2, 2]
+    estr = MEstimator(psi, init=[0., 0.5, 0.5, 0., 0., 0., 0.])
+    estr.estimate(solver='lm')
 
-The ``ee_ipw`` supports both binary and continuous outcomes automatically. Both of these variable types are handled in
-the same way due to the form of the Horwitz-Thompson estimator.
+After successful optimization, we can inspect the estimated values.
+
+.. code::
+
+    estr.theta[0]    # causal mean under X1
+    estr.theta[1:]   # regression coefficients
 
-Unlike the GEE-trick for IPW (which provides a conservative estimator of the variance), the variance estimator here
-is correct. This means it will be narrower than the GEE-trick. Therefore, this approach is generally preferred over
-the GEE-trick to calculating the variance for the IPW estimator. It is also much more computationally efficient than
-the bootstrap.
 
 Augmented inverse probability weighting
 ----------------------------------------------
 
-Before, we model the outcome and treatment models separately. Now, we will consider the augmented inverse probability
-weighting (AIPW) model, which incorporates both the treatment and outcome models. AIPW is a semi-parametric
-doubly-robust estimator for the average treatment effect. For a basic overview, see Funk et al. (2011).
+Finally, we use the augmented inverse probability weighting (AIPW) esitmator, which incorporates both a model for
+:math:`Y` conditional on :math:`A` and :math:`W`, and a model for :math:`A` conditional on :math:`W`. The AIPW estimator
+for the mean difference can be written as
+
+.. math::
+
+    \frac{1}{n} \sum_{i=1}^n \frac{A_i \times Y_i}{\pi_i} - \frac{m_1(W_i; \hat{\beta})(A_i-\pi_i}{\pi_i} -
+    \frac{1}{n} \sum_{i=1}^n \frac{(1-A_i) \times Y_i}{1-\pi_i} + \frac{m_0(W_i; \hat{\beta})(A_i-\pi_i}{1-\pi_i}
 
-*A limitation*: as with g-computation, the built-in AIPW estimating equation only uses a single outcome model
-and that outcome model does *not* support interaction terms. Here the AIPW is meant as a basic example. For
-more general use, the provided estimating equation should be adapted. But the built-in estimating equation will provide
-a basic structure for user's to build off of.
 
-The estimating equations for the AIPW estimator are also provided in ``delicatessen``. To load the estimating equations,
-we call
+where :math:`m_a(W_i; \hat{\beta}) = E[Y_i|A_i=a,W_i; \hat{\beta}]`, and
+:math:`\pi_i = Pr(A_i = 1 | W_i; \hat{\alpha})`. In ``delicatessen``, the built-in AIPW estimator consists of 5
+estimating equations, with both binary and continuous outcomes supported. Similar to IPW (and unlike g-computation),
+the built-in AIPW estimator only supports the average causal effect as the parameter of interest.
+
+The stacked estimating equations are
+
+.. math::
+
+    \sum_i^n \psi_0(Y_i, A_i, \pi_i, \theta_0) = \sum_i^n (\theta_1 - \theta_2) - \theta_0 = 0
+
+    \sum_i^n \psi_1(Y_i, A_i, W_i, \pi_i, \theta_1) = \sum_i^n (\frac{A_i \times Y_i}{\pi_i} -
+    \frac{\hat{Y^1}(A_i-\pi_i}{\pi_i}) - \theta_1 = 0
+
+    \sum_i^n \psi_0(Y_i, A_i, \pi_i, \theta_2) = \sum_i^n (\frac{(1-A_i) \times Y_i}{1-\pi_i} +
+    \frac{\hat{Y^0}(A_i-\pi_i}{1-\pi_i})) - \theta_2 = 0
+
+    \sum_i^n \psi_g(A_i, W_i, \alpha) = \sum_i^n (A_i - expit(W_i^T \alpha)) W_i = 0
+
+    \sum_i^n \psi_m(Y_i, X_i, \beta) = \sum_i^n (Y_i - X_i^T \beta) X_i = 0
+
+where :math:`\theta_1` is the average causal effect, :math:`\theta_2` is the mean under the first plan, :math:`\theta_3`
+is the mean under the second, :math:`\alpha` is the parameters for the propensity score logistic model, and
+:math:`\beta` is the parameters for the regression model used to predict the outcomes. For binary outcomes, the final
+estimating equation is replaced with the logistic model analog.
+
+To load the estimating equations,
 
 .. code::
 
     from delicatessen import MEstimators
     from delicatessen.estimating_equations import ee_aipw
 
-As always, we will wrap the built-in estimating equation inside a function.
+The estimating equation is then wrapped inside the wrapper ``psi`` function. Like ``ee_gformula``, ``ee_aipw`` requires
+some additional data prep. Specifically, we need to create a copy of the data set where :math:`A=1` for everyone and another
+copy where :math:`A=0` for everyone. Below is code that does this step and creates the wrapper function
 
 .. code::
 
-    def psi(theta):
-        return ee_aipw(theta, X=d[['C', 'A', 'W']], y=d['Y'], treat_index=1)
+    # Creating data under the plans
+    d1 = d.copy()
+    d1['A'] = 1
+    d0 = d.copy()
+    d0['A'] = 0
 
-The arguments for ``ee_aipw`` are the :math:`\theta` values, the covariates (including an intercept (C) and the
-action (A)), the outcome values (Y), and the column index for the treatment in X. Here, 1 designates the second column
-(python uses zero-indexing), which corresponds to ``'A'`` in how the X data is formatted.
+    # Creating interaction terms
+    d['AxW'] = d['A'] * d['W']
+    d1['AxW'] = d1['A'] * d1['W']
+    d0['AxW'] = d0['A'] * d0['W']
 
-Now we can call the M-estimator to solve for the values and the variance. Here, the initial values provided must be
-3+*b*+*b-1* (where *b* is the number of columns in X). This is because the AIPW estimating equations output the
-average treatment effect, risk under all-treated, risk under none-treated, and the outcome model coefficients, and
-the treatment model coefficients.
+    def psi(theta):
+        return ee_gformula(theta,                        # Parameters
+                           y=d['Y'],                     # Outcome
+                           A=d['A'],                     # Action
+                           W=d[['C', 'W']],              # PS model
+                           X=d[['C', 'A', 'W', 'AxW']],  # Outcome model
+                           X=d1[['C', 'A', 'W', 'AxW']], # Plan A=1
+                           X=d0[['C', 'A', 'W', 'AxW']]) # Plan A=0
 
-As for starting values, it will likely be best practice to have the initial values set as ``[0., 0.5, 0.5, ...]`` in
-general. The regression initial values can also be pre-washed to speed up optimization.
+Note that we add an intercept to the outcome model by adding a column of 1's via ``d['C']``.
 
-.. code::
+Here, the initial values provided must be 3+*b*+*c* (where *b* is the number of columns in W and *c* is the number of
+columns in X). For binary outcomes, it will likely be best practice to have the initial values set as
+``[0., 0.5, 0.5, ...]``. followed by b ``0.``'s. For continuous outcomes, all ``0.`` can be used instead. Furthermore,
+a regression models could be optimized outside of ``delicatessen`` and those (pre-washed) regression estimates can be
+passed as initial values to speed up optimization.
 
-    mestimation = MEstimator(stacked_equations=psi, init=[0., 0.5, 0.5, 0., 0., 0.])
-    mestimation.estimate(solver='lm')
-
-Now the average treatment effect, as well as the variance, can be output. Here, a key advantage of M-estimation can be
-seen. The form of an M-estimator allows us to estimate the variance directly, while appropriately allowing for the
-uncertainty in the regression model parameters to be carried forward. M-estimation does this automatically for us.
-Essentially, we do not need to bootstrap or use the GEE-trick for IPW to estimate the variance!
+Now we can call the M-estimator to solve for the values and the variance.
 
 .. code::
 
-    mestimation.theta[0]
-    mestimation.variance[0, 0]
+    estr = MEstimator(psi, init=[0., 0.5, 0.5,
+                                 0., 0.,
+                                 0., 0., 0., 0.])
+    estr.estimate(solver='lm')
 
-Besides the average treatment effect, the risk / mean under all-treated can be extracted by
+After successful optimization, we can inspect the estimated values.
 
 .. code::
 
-    mestimation.theta[1]
-    mestimation.variance[1, 1]
+    estr.theta[0]     # causal mean difference of 1 versus 0
+    estr.theta[1]     # causal mean under A=1
+    estr.theta[2]     # causal mean under A=0
+    estr.theta[3:5]   # propensity score regression coefficients
+    estr.theta[5:]    # outcome regression coefficients
 
-and the risk / mean under none-treated can be extracted by
+The variance and Wald-type confidence intervals can also be output via
 
 .. code::
 
-    mestimation.theta[2]
-    mestimation.variance[2, 2]
+    estr.variance
+    estr.confidence_intervals()
 
-The variance estimator in this case will match the influence function estimator of the variance that is commonly used
-for AIPW. See Boos & Stefanski (2013) for more detailed discussion on the relation between M-estimation and influence
-curves.
+Here, the M-Estimation sandwich variance is the same as the influence-curve-based variance estimator. Either of these
+approaches correctly incorporates the uncertainty in estimation of the outcome model into the parameter(s) of interest
+(e.g., average causal effect). Therefore, we do not have to rely on the nonparametric bootstrap (computationally
+cumbersome).
 
-Further Readings
-=============================
+
+References and Further Readings
+===============================
 Boos DD, & Stefanski LA. (2013). M-estimation (estimating equations). In Essential Statistical Inference
 (pp. 297-337). Springer, New York, NY.
 
+Cole SR, & Hernán MA. (2008). Constructing inverse probability weights for marginal structural models.
+*American Journal of Epidemiology*, 168(6), 656-664.
+
 Funk MJ, Westreich D, Wiesen C, Stürmer T, Brookhart MA, & Davidian M. (2011). Doubly robust estimation of causal
 effects. *American Journal of Epidemiology*, 173(7), 761-767.
 
+Hernán MA, & Robins JM. (2006). Estimating causal effects from epidemiological data.
+*Journal of Epidemiology & Community Health*, 60(7), 578-586.
+
 Huber PJ. (1992). Robust estimation of a location parameter. In Breakthroughs in statistics (pp. 492-518).
 Springer, New York, NY.
 
-Satten GA, & Datta S. (2001). The Kaplan–Meier estimator as an inverse-probability-of-censoring weighted average.
-*The American Statistician*, 55(3), 207-210.
+Inderjit, Streibig JC & Olofsdotter M. (2002). Joint action of phenolic acid mixtures and its significance in
+allelopathy research. *Physiol Plant* 114, 422–428.
+
+Snowden JM, Rose S, & Mortimer KM. (2011). Implementation of G-computation on a simulated data set: demonstration
+of a causal inference technique. *American Journal of Epidemiology*, 173(7), 731-738.
diff --git a/docs/Custom Equations.rst b/docs/Custom Equations.rst
index adb1d6d..1f9bf67 100644
--- a/docs/Custom Equations.rst
+++ b/docs/Custom Equations.rst
@@ -6,11 +6,11 @@ be specified. Basically, it will allow for any estimating equation to be passed
 equation(s) *must* be unbiased for the theory behind M-estimation to hold). Here, I provide an overview and tips for
 how to build your own estimating equation.
 
-In general, it will be best if you find an explicit paper or book (most likely written by a statistician) that directly
-provides the estimating equation(s) to you. Deriving your own *unbiased* estimating equation may be a lot of effort
-and will require a statistical proof. This section does *not* address this part of M-estimation. Rather, this section
-provides information on how to construct an estimating equation within ``delicatessen``. ``delicatessen`` assumes you
-are giving it a valid estimating equation.
+In general, it will be best if you find an explicit paper or book that directly provides the estimating equation(s) to
+you. Alternatively, if you can find the score function or gradient for a regression model, that is the corresponding
+estimating equation. This section does *not* address how to derive your own  estimating equation(s). Rather, this
+section provides information on how to construct an estimating equation within ``delicatessen``, as ``delicatessen``
+assumes you are giving it a valid estimating equation.
 
 Building from scratch
 -------------------------------------
@@ -19,7 +19,7 @@ First, we will go through the case of building an estimating equation completely
 go through an example with linear regression. This is how I went about building the ``ee_linear_regression``
 functionality.
 
-First, we have the estimating equation (which is the score function in this case) provided in Boos & Stefanski (2013)
+First, we have the estimating equation (which is the score function) provided in Boos & Stefanski (2013)
 
 .. math::
 
@@ -39,7 +39,8 @@ We will demonstrate using the following simulated data set
 
 
 First, we can build the estimating equation using a for-loop where each ``i`` piece will be stacked together. While this
-for-loop approach will be slow, it is often a good strategy to implement this version first.
+for-loop approach will be slow, it is often a good strategy to implement a for-loop version that is easier to debug
+first.
 
 Below calculates the estimating equation for each ``i`` in the for-loop. This function returns a stacked array of each
 ``i`` observation as a 3-by-n array. That array can be validly passed to the ``MEstimator`` for optimization and
@@ -73,13 +74,13 @@ We can then run this estimating equation with
     mest = MEstimator(psi, init=[0., 0., 0.])
     mest.estimate()
 
-for which the coefficients match the coefficients from a ordinary least squares model (variance estimates will differ).
-Here, we can further vectorize the estimating equation. In the vector-form, this code will run much faster and this
-is often the best approach to boosting speed in terms of run-time.
+for which the coefficients match the coefficients from a ordinary least squares model (variance estimates will differ,
+since most OLS software uses a different variance estimator). Here, we can further vectorize the estimating equation. In
+the vector-form, this code will run much faster.
 
 With some careful experimentation, the following is a vectorized version. Remember that ``delicatessen`` is expecting a
-b-by-n array to be given by the ``psi`` function. Failure to provide this is a common mistake when building custom
-estimating equations.
+3-by-n array to be given by the ``psi`` function in this example. Failure to provide this is a common mistake when
+building custom estimating equations.
 
 .. code::
 
@@ -92,8 +93,8 @@ estimating equations.
 
 As before, we can run this chunk of code. However, this is substantially faster. If we run both implementations on the
 same data set of 10,000 observations, the for-loop version took approximately 1.60 seconds and the vectorized version
-took 0.05 seconds (on my fairly new laptop). That is a large difference in run time! Vectorizing (even parts of an
-estimating equation) can help to improve run-times if you find the M-Estimation procedure taking too long.
+took 0.05 seconds (on my fairly new laptop). Vectorizing (even parts of an estimating equation) can help to improve
+run-times if you find the M-Estimation procedure taking too long.
 
 
 Building with basics
@@ -102,19 +103,19 @@ Building with basics
 Instead of building everything from scratch, you can also piece together the built-in estimating equations with your
 own code. To demonstrate this, I will go through how I developed the code for inverse probability weighting.
 
-The inverse probability weighting estimator consists of four estimating equations: the propensity score model, the
-weighted mean for treatment, the weighted mean for no-treatment, and the difference between the weighted means. We
+The inverse probability weighting estimator consists of four estimating equations: the difference between the weighted
+means, the weighted mean under :math:`A=1`, the weighted mean under :math:`A=0`, and the propensity score model. We
 can write this as
 
 .. math::
 
-    \sum_i^n \psi_1(A_i, W_i, \theta) = \sum_i^n (A_i - expit(W_i^T \theta)) W_i = 0
+    \sum_i^n \psi_1(Y_i, A_i, \pi_i, \theta_0) = \sum_i^n (\theta_1 - \theta_2) - \theta_0 = 0
 
     \sum_i^n \psi_2(Y_i, A_i, \pi_i, \theta_1) = \sum_i^n \frac{A_i \times Y_i}{\pi_i} - \theta_1 = 0
 
     \sum_i^n \psi_3(Y_i, A_i, \pi_i, \theta_2) = \sum_i^n \frac{(1-A_i) \times Y_i}{1-\pi_i} - \theta_2 = 0
 
-    \sum_i^n \psi_4(Y_i, A_i, \pi_i, \theta_0) = \sum_i^n (\theta_1 - \theta_2) - \theta_0 = 0
+    \sum_i^n \psi_4(A_i, W_i, \alpha) = \sum_i^n (A_i - \text{expit}(W_i^T \alpha)) W_i = 0
 
 
 Rather than re-code the logistic regression model (to estimate the propensity scores), we will use the built-in
@@ -130,13 +131,13 @@ estimator
         y = np.asarray(y)
         beta = theta[3:]   # Extracting out theta's for the regression model
 
-        # Estimating propensity score
+        # Estimating propensity score using delicatessen
         preds_reg = ee_logistic_regression(theta=beta,    # Using logistic regression
                                            X=W,           # Plug-in covariates for X
                                            y=A)           # Plug-in treatment for Y
 
         # Estimating weights
-        pi = inverse_logit(np.dot(W, beta))          # Getting Pr(A|W) from model (using delicatessen.utilities)
+        pi = inverse_logit(np.dot(W, beta))          # Pr(A|W) using delicatessen.utilities
 
         # Calculating Y(a=1)
         ya1 = (A * y) / pi - theta[1]                # i's contribution is (AY) / \pi
@@ -176,5 +177,6 @@ Here is a list of common mistakes, most of which I have done myself.
 4. The ``theta`` values and ``b`` *must* be in the same order. If ``theta[0]`` is the mean, the 1st row of the returned
    array better be the mean!
 
-If you still have trouble, please open an issue on `GitHub<https://github.com/pzivich/Delicatessen/issues>`_. This will
-help me to add other common mistakes here and improve the documentation for custom estimating equations.
+If you still have trouble, please open an issue at
+`pzivich/Delicatessen <https://github.com/pzivich/Delicatessen/issues>`_. This will help me to add other common
+mistakes here and improve the documentation for custom estimating equations.
diff --git a/docs/Examples.rst b/docs/Examples.rst
index 540ff63..0bf7ae9 100644
--- a/docs/Examples.rst
+++ b/docs/Examples.rst
@@ -1,9 +1,11 @@
 Examples
 =====================================
 
-Here, I will go through a variety of examples (using custom estimating equations, rather than the built-in options).
-These examples follow Chapter 7 of Boos & Stefanski. If you have the book (or access to it), then reading along with
-each section may be helpful. To make it easier, I also provide the section heading from Chapter 7.
+Here, we will implement some of the examples described in Chapter 7 of Boos & Stefanski (2013). If you have the book
+(or access to it), then reading along with each section may be helpful. To make it easier, I also provide the section
+heading from Chapter 7. Finally, we will code each of the estimating equations by-hand (rather than using the built-in
+options).
+
 
 Boos DD, & Stefanski LA. (2013). M-estimation (estimating equations). In Essential Statistical Inference
 (pp. 297-337). Springer, New York, NY.
@@ -12,22 +14,34 @@ Boos DD, & Stefanski LA. (2013). M-estimation (estimating equations). In Essenti
 Sample Mean (7.2.2)
 -------------------------------
 
-First, we demonstrate a simple estimating equation for the mean and variance (stacked together). The estimating
-equations are
+First, we demonstrate a simple estimating equation for the mean and variance. The estimating equations are
 
 .. math::
 
     \psi_1(Y_i, \theta) &= Y_i - \theta_1 \\
     \psi_2(Y_i, \theta) &= (Y_i - \theta_1)^2 - \theta_2
 
-To demonstrate the example, we will use some generic data for Y
+To demonstrate the example, we will use some generic data for :\math:`Y`. Below is an example data set that will be
+used up to Section 7.2.6:
 
 .. code::
 
     import numpy as np
-    y = np.array([-2, 1, 3, 4, 1, 4, -5, 3, 6])
 
-For use with the M-estimator in ``delicatessen``, we can program these estimating equations as
+    np.random.seed(80950841)
+
+    n = 200
+    data = pd.DataFrame()
+    data['Y'] = np.random.normal(loc=10, scale=2, size=n)
+    data['X'] = np.random.normal(loc=5, size=n)
+    data['C'] = 1
+
+    # Subsetting the data for 7.2.2
+    y = np.asarray(data['Y'])
+
+
+For use with the M-estimator in ``delicatessen``, we need to write the Python analog of the estimating equation math.
+Below is one way we can do this. Notice that the shape of the returned array in 2-by-n.
 
 .. code::
 
@@ -36,53 +50,62 @@ For use with the M-estimator in ``delicatessen``, we can program these estimatin
         vari = (y-theta[0])**2 - theta[1]   # \psi_2 from above
         return mean, vari
 
-After defining the stacked estimating equations, the M-estimator can be called
+Once we have written up our stacked estimating equations, ``MEstimator`` can be called to solve for :math:`\theta` and
+the sandwich variance. We can do that via
 
 .. code::
 
-    mestimate = MEstimator(psi, init=[0, 1])
-    mestimate.estimate()
+    from delicatessen import MEstimator
+
+    estr = MEstimator(psi, init=[0, 1])
+    estr.estimate()
+
+    print(estr.theta)                # [10.163  4.112]
+    print(estr.asymptotic_variance)  # [[ 4.112, -1.674], [-1.674, 36.164]]
+
 
-The M-estimator will solve for :math:`\theta` via a root finding procedure. For the sandwich variance, ``delicatessen``
-uses a numerical approximation procedure for the derivative. This is different from the closed-form variance estimator
-provided in Chapter 7, but both should return the same answer (within computational error). The advantage of the
-numerical derivatives are that they can be done for arbitrary estimating equations.
+The M-Estimator solves for :math:`\theta` via a root finding procedure given the initial values in ``init``. Since the
+variance must be >0, we provide a positive initial value. For the sandwich variance, ``delicatessen`` uses a numerical
+approximation procedure for the derivative. This is different from the closed-form variance estimator provided in
+Chapter 7, but both should return the same answer (within computational error). The advantage of the numerical
+derivatives are that they can be done for arbitrary estimating equations.
 
 
 Ratio (7.2.3)
 -------------------------------
 
-For the ratio estimator, the following estimating equation is provided
+Now consider if we wanted to estimate the ratio between two means. For estimation of a ratio, we can consider the
+following estimating equation
 
 .. math::
 
-    \psi_1(Y_i, X_i, \theta) = Y_i - \theta_1 \X_i
+    \psi_1(Y_i, X_i, \theta) = Y_i - \theta_1 X_i
 
-For use with the M-estimator in ``delicatessen``, we can program this estimating equation as
 
-.. code::
-
-    import numpy as np
-    import pandas as pd
+We can translate the estimating equation from math into python as
 
-    n = 200
-    data = pd.DataFrame()
-    data['Y'] = np.random.normal(loc=10, scale=2, size=n)
-    data['X'] = np.random.normal(loc=5, size=n)
+.. code::
 
     def psi(theta):
-        return data['Y'] - data['X']*theta
+        return data['Y'] - data['X']*theta  # \psi_1 from above
+
 
-As before, we can solve the estimating equation by
+Now, we can pass this estimating equation and data to `MEstimator`
 
 .. code::
 
-    mestimate = MEstimator(psi, init=[1, ])
-    mestimate.estimate()
+    estr = MEstimator(psi, init=[1, ])
+    estr.estimate()
 
-As you may notice, only a single initial value is provided (since only a single array is being returned).
+    print(estr.theta)                # [2.082]
+    print(estr.asymptotic_variance)  # 0.338
 
-The chapter also provides an alternative series of estimating equations for the ratio.
+As you may notice, only a single initial value is provided (since only a single array is being returned). Furthermore,
+we provide an initial value >0 since we are estimating a ratio.
+
+There is another set of stacked estimating equations we can consider for the ratio. Specifically, we can estimate each
+of the means and then take the ratio of those means (rather than doing everything simultaneously). Below is this
+alternative set of estimating equations
 
 .. math::
 
@@ -91,58 +114,69 @@ The chapter also provides an alternative series of estimating equations for the
     \psi_3(Y_i, X_i, \theta) &= \theta_1 - \theta_2 \times \theta_3
 
 
-Similarly, this can also be programmed via
+Translating this to an estimating equation in Python
 
 .. code::
 
     def psi(theta):
-        mean_y = data['Y'] - theta[0]
-        mean_x = data['X'] - theta[1]
-        ratio = np.ones(data.shape[0]) * (theta[0] - theta[1]*theta[2])
+        mean_y = data['Y'] - theta[0]        # \psi_1 from above
+        mean_x = data['X'] - theta[1]        # \psi_2 from above
+        ratio = (np.ones(data.shape[0]) *    # \psi_3 from above
+                 (theta[0] - theta[1]*theta[2]))
         return mean_y, mean_x, ratio
 
-    mestimate = MEstimator(psi, init=[0, 0, 1])
-    mestimate.estimate()
+    estr = MEstimator(psi, init=[0, 0, 1])
+    estr.estimate()
+
+    print(estr.theta)  # [10.163,  4.880,  2.082]
 
-Here, there is a series of estimating equations. It is also important to note the use of ``np.ones`` in the third step.
-This ensures that ``ratio`` consists on *n* observations. Without multiplying by the array of ones, ``ratio`` would be a
-single value. However, ``MEstimator`` expects a 3-by-*n* array here. Multiplying the 3rd equation by an array of 1's
-keeps the correct dimension and keeps the values.
+Here, we used a trick to make sure the dimension of ``ratio`` stays as :math:`n`, we use ``np.ones``. Without
+multiplying by the array of ones, ``ratio`` would be a single value. However, ``MEstimator`` expects a
+:math:`3 \times n` array. Multiplying the 3rd equation by an array of 1's keeps the correct dimension.
 
-Also notice that these estimating equations require the use of 3 ``init`` values, unlike the other ratio estimator.
+Also notice this form requires the use of 3 ``init`` values, unlike the other ratio estimator. As before, the ratio
+initial value is set >0 to be nice to the root-finder.
 
 
 Delta Method (7.2.4)
 -------------------------------
 
-The delta method can also be cast as an M-estimation problem. The chapter demonstrates two transformations of Y and
-their corresponding mean. The stacked estimating equations are
+The delta method has been used in a variety of contexts, including estimating the variance for transformations of
+parameters. Instead of separately estimating the parameters, transforming the parameters, and then using the delta
+method to estiamte the variance of the transformed parameters; we can apply the transformation in an estimating
+equation and automatically estimate the variance for the transformed parameter(s) via the sandwich variance. To do this,
+we stack the estimating equation for the transformation into our set of estimating equations. Below is the
+mean-variance estimating equations stacked with two transformations of the variance
 
 .. math::
 
     \psi_1(Y_i, \theta) &= Y_i - \theta_1 \\
     \psi_2(Y_i, \theta) &= (Y_i - \theta_1)^2 - \theta_2 \\
     \psi_3(Y_i, \theta) &= \sqrt{\theta_2} - \theta_3 \\
-    \psi_4(Y_i, \theta) &= \log(\theta_2) - \theta_4 \\
+    \psi_4(Y_i, \theta) &= \log(\theta_2) - \theta_4
 
 
-These equations can be expressed programmatically for ``delicatessen`` as
+These equations can be expressed programmatically as
 
 .. code::
 
     def psi_delta(theta):
-        return (data['Y'] - theta[0],
-                (data['Y'] - theta[0])**2 - theta[1],
-                np.ones(data.shape[0])*np.sqrt(theta[1]) - theta[2],
-                np.ones(data.shape[0])*np.log(theta[1]) - theta[3])
+        mean = data['Y'] - theta[0]                           # \psi_1 from above
+        variance = (data['Y'] - theta[0])**2 - theta[1]       # \psi_2 from above
+        sqrt_var = (np.ones(data.shape[0])*np.sqrt(theta[1])  # \psi_3 from above
+                    - theta[2])
+        log_var = (np.ones(data.shape[0])*np.log(theta[1])    # \psi_4 from above
+                   - theta[3])
+        return (mean, variance, sqrt_var, log_var)
 
-Notice the use of the ``np.ones`` trick as done with the ratio estimating equations to ensure that the final equations are
-the correct shapes.
+Notice the use of the ``np.ones`` trick again to ensure that the final equations are the correct shapes.
 
 .. code::
 
-    mestimate = MEstimator(psi, init=[0, 0, 1, 1])
-    mestimate.estimate()
+    estr = MEstimator(psi, init=[0, 1, 1, 1])
+    estr.estimate()
+
+    print(estr.theta)  # [10.163, 4.112, 2.028, 1.414]
 
 Here, there are 4 stacked equations, so ``init`` must be provided 4 values.
 
@@ -150,25 +184,31 @@ Here, there are 4 stacked equations, so ``init`` must be provided 4 values.
 Instrumental Variable (7.2.6)
 -------------------------------
 
-The first set of estimating equations for the instrumental variable analysis are
+As a further example, consider the following instrumental variable approach to correcting for measurement error of a
+variable. Here, :math:`Y` is the outcome of interest, :math:`X` is the true, unmeasured variable, :math:`W` is the
+possibly mismeasured variables, and :math:`T` is the instrument for :math:`X`.
+
+The first set of estimating equations consider in Chapter 7 are
 
 .. math::
 
-    \psi_1(Y_i, W_i, T_i, \theta) &= \theta_1 - T \\
-    \psi_2(Y_i, W_i, T_i, \theta) &= (Y - \theta_2 W)(\theta_1 - T) \\
+    \psi_1(Y_i, W_i, T_i, \theta) &= \theta_1 - T_i \\
+    \psi_2(Y_i, W_i, T_i, \theta) &= (Y_i - \theta_2 W_i)(\theta_1 - T_i)
 
-To demonstrate the example, below is some generic simulated data
+To demonstrate the example, below is some generic simulated data in the described instrumental variable context
 
 .. code::
 
+    np.random.seed(809421)
     n = 500
+
     data = pd.DataFrame()
     data['X'] = np.random.normal(size=n)
     data['Y'] = 0.5 + 2*data['X'] + np.random.normal(loc=0, size=n)
     data['W'] = data['X'] + np.random.normal(loc=0, size=n)
     data['T'] = -0.75 - 1*data['X'] + np.random.normal(loc=0, size=n)
 
-These estimating equations can be programmed for ``delicatessen`` as
+The previous estimating equations can be translated as
 
 .. code::
 
@@ -176,20 +216,22 @@ These estimating equations can be programmed for ``delicatessen`` as
         return (theta[0] - data['T'],
                 (data['Y'] - data['W']*theta[1])*(theta[0] - data['T']))
 
-    mestimate = MEstimator(psi, init=[0.1, 0.1, ])
-    mestimate.estimate()
+    estr = MEstimator(psi, init=[0.1, 0.1, ])
+    estr.estimate()
+
+    print(estr.theta)  # [-0.777,  1.769,]
 
-As mentioned in the chapter, certain joint distributions may be of interest. To capture those distributions, the
-estimating equations from before were further updated to
+As mentioned in the chapter, certain joint distributions may be of interest. To capture these additional distributions,
+the estimating equations were updated to
 
 .. math::
 
-    \psi_1(Y_i, W_i, T_i, \theta) &= \theta_1 - T \\
-    \psi_2(Y_i, W_i, T_i, \theta) &= \theta_2 - W \\
-    \psi_3(Y_i, W_i, T_i, \theta) &= (Y - \theta_3 W)(\theta_2 - W) \\
-    \psi_4(Y_i, W_i, T_i, \theta) &= (Y - \theta_4 W)(\theta_1 - T) \\
+    \psi_1(Y_i, W_i, T_i, \theta) &= \theta_1 - T_i \\
+    \psi_2(Y_i, W_i, T_i, \theta) &= \theta_2 - W_i \\
+    \psi_3(Y_i, W_i, T_i, \theta) &= (Y_i - \theta_3 W_i)(\theta_2 - W_i) \\
+    \psi_4(Y_i, W_i, T_i, \theta) &= (Y_i - \theta_4 W_i)(\theta_1 - T_i)
 
-Again, we can easily write these equations for ``delicatessen``,
+Again, we can easily translate these equations for ``delicatessen``,
 
 .. code::
 
@@ -200,40 +242,151 @@ Again, we can easily write these equations for ``delicatessen``,
                 (data['Y'] - data['W']*theta[3])*(theta[0] - data['T'])
                 )
 
-    mestimator = MEstimator(psi, init=[0.1, 0.1, 0.1, 0.1])
-    mestimator.estimate()
+    estr = MEstimator(psi, init=[0.1, 0.1, 0.1, 0.1])
+    estr.estimate()
 
-This example further demonstrates the flexibility of M-estimation by stacking together estimating equations.
+    print(estr.theta)  # [-0.777, 0.005, 0.964, 1.769]
+
+This example demonstrates the flexbility of M-Estimation through the ability to stack estimating equations together.
 
 
 Robust Location (7.4.1)
 -------------------------------
 
-For robust location estimation, the estimating equation is
+To begin, we generate some generic data used for this example and several of the following
+
+.. code::
+
+    np.random.seed(7841)
+    y = np.random.normal(size=250)
+    n = y.shape[0]
+
+For the robust mean, the estimating equation proposed by Huber (1964) is
+
+.. math::
+
+    \psi_k(Y_i, \theta) = g_k(Y_i) - \theta_1
+
+where :math:`k` indicates the bound, such that if :math:`Y_i>k` then :math:`k`, or :math:`Y_i<-k` then :math:`-k`,
+otherwise :math:`Y_i`. Below is the estimating equation translated into code
+
+.. code::
+
+    def psi_robust_mean(theta):
+        k = 3                          # Bound value
+        yr = np.where(y > k, k, y)     # Applying upper bound
+        yr = np.where(y < -k, -k, y)   # Applying lower bound
+        return yr - theta
+
+
+    estr = MEstimator(psi_robust_mean, init=[0.])
+    estr.estimate()
+
+    print(estr.theta)  # [-0.0126]
+
+Notice that the estimating equation here is not smooth. Specifically, there is a jump at :math:`k`. Therefore, this
+estimator only behaves correctly when the values of :math:`\theta` are differentiable (i.e., the true mean can't be
+at :math:`k` or :math:`-k`).
+
+
+Quantile Estimation (7.4.2)
+-------------------------------
+
+Despite the sandwich variance needing the function to be smooth at :math:`\theta` (so it is differentiable),
+M-Estimation can also be used with non-smooth function. For example, the estimating equations for the sample quantile
+is
 
 .. math::
 
-    \psi_k(Y_i, \theta) = Y^k_i - theta_1
+    \psi_q(Y_i, \theta) = q - I(Y_i \le \theta)
+
+It is this section, that we need to talk about different root-finding methods, and numerically approximating
+derivatives. In the previous examples, we had smooth function that were both easy to find the roots of and had smooth
+functions for derivatives. However, that is not the case for quantile estimation. So, we need to use some 'tricks' to
+help the procedure along.
 
-where *k* indicates the upper and lower bound, and Y superscript *k* is the bounded values of Y.
+First, we are going to use the ``'hybr'`` method. we have found this method to be more reliable when attempting to find
+the roots. Often the ``'lm'`` and ``'newton'`` methods appear worse at exploring the space. Next, our estimating
+equations 'jump' in terms of their returned value (i.e., they are not smooth). This comes in to the ``tolerance``
+parameter. The tolerance determines whether the root-finding has converged. For many quantiles were aren't going to
+reach the strict tolerance values. So, we are going to weaken them (the algorithm will be considered as converged under
+a weaker condition). If this is not changed, then a non-convergence error will come back.
 
-Below is the estimating equation in Python
+Now we can talk about numerically approximating the derivatives. Numerical approximations roughly work by calculating
+the slope of a line from two points on either side of value (akin to the definition of a derivative you may remember
+from math class). For smooth functions, we can choose these points 'close' to the true value. However, this is not the
+case for non-smooth functions. For non-smooth functions the derivative can be poorly approximated when relying on points
+'too close' to the value. We can address this issue by increasing the ``dx`` parameter. However, large ``dx`` parameters
+can also lead to poor approximations. Therefore, we will also increase the ``order`` parameter, which controls the
+number of points to use (note: it must be odd).
+
+Now, that we have these tricks, we are ready to find the 25th, 50th, and 75th percentiles using M-Estimation. The
+estimating equations are
 
 .. code::
 
-    import numpy as np
-    var = np.array([1, -10, 2, 1, 4, 1, 4, 2, 4, 2, 3, 12])
+    def psi_quantile(theta):
+        return (0.25 - 1*(y <= theta[0]),
+                0.50 - 1*(y <= theta[1]),
+                0.75 - 1*(y <= theta[2]),)
 
-    def psi(theta):
-        var = np.where(var > k, k, var)       # Apply the upper bound
-        var = np.where(var < -k, -k, var)     # Apply the lower bound
-        return var - theta                    # Estimating equation for robust mean
 
-    mestimator = MEstimator(psi, init=[0., ])
-    mestimator.estimate()
+    estr = MEstimator(psi_quantile, init=[0., 0., 0.])
+    estr.estimate(solver='hybr',   # Selecting the hybr method
+                  tolerance=1e-3,  # Increasing the tolerance
+                  dx=1,            # Increasing distance for numerical approx
+                  order=9)         # Increasing the number of points for numerical approx
+
+    print(estr.theta)  # [-0.597  0.048  0.740]
+
+We can compare these values to
+
+.. code::
+
+    np.quantile(y, q=[0.25, 0.50, 0.75])  # [-0.592, 0.047, 0.740]
+
+You'll notice that there is a slight difference. This difference is a result of the non-smooth function. Values 'close'
+to these points will not improve the zero finding in the estimating equations. That was why we decreased the tolerance
+originally. So, there may be a slight discrepancy between the closed-form solution and M-Estimation.
+
+For non-smooth functions, it is good practice to check against some closed form for the estimating equations.
+
+
+Positive Mean Deviation (7.4.3)
+-------------------------------
+
+For another non-smooth estimating equation(s), we can talk about the positive mean deviation. The estimating equations
+are
+
+.. math::
+
+    \psi_1 = 2(Y_i - \theta_2) I(Y_i > \theta_2) - \theta_1 \\
+    \psi_2 = 0.5 - I(Y_i \le \theta_2)
+
+where :math:`\theta_1` is the positive mean deviation and :math:`\theta_2` is the median.
+
+The estimating equations can be translated into code by
+
+.. code::
+
+    def psi_deviation(theta):
+        return ((2 * (y - theta[1]) * (y > theta[1])) - theta[0],
+                1/2 - (y <= theta[1]), )
+
+As before, we will use the ``'hybr'`` method along with the updated parameters
+
+.. code::
+
+    estr = MEstimator(psi_deviation, init=[0., 0., ])
+    estr.estimate(solver='hybr',   # Selecting the hybr method
+                  tolerance=1e-3,  # Increasing the tolerance
+                  dx=1,            # Increasing distance for numerical approx
+                  order=9)         # Increasing the number of points for numerical approx
+
+    print(estr.theta)  # [0.803 0.042]
 
-Notice that the estimating equation here is not smooth. Specifically, there is a jump at *k*. Therefore, this will only
-work for values of theta that are differentiable (i.e., the true mean can't be at *k*).
+If we had used the closed-form definition, we would have ended up with (0.798, 0.047). These values are close, and
+again due to the non-smooth nature of the estimating equations.
 
 
 Linear Regression (7.5.1)
@@ -243,16 +396,14 @@ For linear regression, the estimating equation is
 
 .. math::
 
-    \psi(X_i, Y_i \beta) = (Y_i - X_i^T \beta) X_i
+    \psi(X_i, Y_i, \beta) = (Y_i - X_i^T \beta) X_i
 
-Here, we present the vectorized version first. Notice that the vectorized version will generally be faster than a
-for-loop implementation. However, it is easy to make a mistake with a vectorized version, so I generally recommend
-creating a for-loop version first (and then creating a vectorized version if that for-loop is too slow).
 
-With some generic data,
+For the following examples, the following generic simulated data is used
 
 .. code::
 
+    np.random.seed(5555)
     n = 500
     data = pd.DataFrame()
     data['X'] = np.random.normal(size=n)
@@ -260,7 +411,11 @@ With some generic data,
     data['Y'] = 0.5 + 2*data['X'] - 1*data['Z'] + np.random.normal(loc=0, size=n)
     data['C'] = 1
 
-The estimating equation and M-estimation procedure is then called by
+As with all the preceding estimating equations, there are multiple ways to code these. Since linear regression involes
+some careful matrix manipulations for the programmed estimating equations to return the correct format for
+``delicatessen``, we highlight two variations here.
+
+First, we present a vectorized version first.
 
 .. code::
 
@@ -270,10 +425,13 @@ The estimating equation and M-estimation procedure is then called by
         beta = np.asarray(theta)[:, None]
         return ((y - np.dot(x, beta)) * x).T
 
-    mestimator = MEstimator(psi_regression, init=[0.1, 0.1, 0.1])
-    mestimator.estimate()
+    estr = MEstimator(psi_regression, init=[0., 0., 0.])
+    estr.estimate()
+
+    print(estr.theta)  # [0.477, 2.123, -0.852]
 
-As mentioned, a for-loop version can also be used. Below is an example of the for-loop version for regression
+For the second approach, a for-loop variation is used instead. Below is the for-loop equivalent for the estimating
+equations
 
 .. code::
 
@@ -296,8 +454,52 @@ As mentioned, a for-loop version can also be used. Below is an example of the fo
         return np.asarray(est_vals).T
 
 
-GEE (7.5.6)
+While these two approaches give the same answer, vectorized versions will generally be faster than for-loop variations
+(but may be less 'human readable'). For example, the vectorized version has a runtime of 0.037 seconds and the for-loop
+version has a runtime of 0.166 seconds (4.5x slower) on my laptop. Having said that, it is easy to make a mistake with
+a vectorized version. I would generally recommend creating a for-loop version first (and then creating a vectorized
+version if that for-loop is too slow for your purposes).
+
+
+Robust Regression (7.5.4)
 -------------------------------
 
-... to be added ...
+The next example is robust regression, where the standard linear regression model is made robust to outliers.
+Essentially, we use the robust mean formula from before but now apply it to the error terms of the regression model.
+The estimating equations are
+
+.. math::
+
+    \psi(X_i, Y_i \beta) = g_k(Y_i - X_i^T \beta) X_i
+
+where :math:`k` indicates the bound, such that if :math:`Y_i>k` then :math:`k`, or :math:`Y_i<-k` then :math:`-k`,
+otherwise :math:`Y_i`.
+
+Taking the previous vectorized version of the linear regression model and building in the :math:`g_k()` function,
+
+.. code::
+
+    def psi_regression(theta):
+        X = np.asarray(data[['C', 'X', 'Z']])
+        y = np.asarray(data['Y'])[:, None]
+        beta = np.asarray(theta)[:, None]
+        k = 2
+
+        # Generating predictions and applying Huber function for robust
+        preds = np.asarray(y - np.dot(X, beta))
+        preds = np.where(preds > k, k, preds)       # Apply the upper bound
+        preds = np.where(preds < -k, -k, preds)     # Apply the lower bound
+
+        # Output b-by-n matrix
+        return (preds * X).T
+
+
+    estr = MEstimator(psi_regression, init=[0., 0., 0.])
+    estr.estimate()
+
+    print(estr.theta)  # [0.491, 2.05, -0.795]
+
+You'll notice that the coefficients have changed slightly here. That is because we have reduced the extent of outliers
+on the estimation of the linear regression parameters (however, our simulated data mechanism doesn't really result in
+major outliers, so the change is small here).
 
diff --git a/docs/Optimization Advice.rst b/docs/Optimization Advice.rst
index af2fbd1..1a5c4f7 100644
--- a/docs/Optimization Advice.rst
+++ b/docs/Optimization Advice.rst
@@ -55,11 +55,18 @@ By default, ``delicatessen`` uses the secant method available in ``scipy.optimiz
 also supports algorithms in ``scipy.optimize.root``. I personally have had the most success with the 'secant' in, but
 I have also seen the 'lm' optimization method perform well (some internal tests check with lm as well).
 
+Non-smooth equations
+--------------------
+As mentioned in the examples, non-smooth estimating equations (e.g., percentiles, positive mean deviation, etc.) can be
+difficult to optimize. In these cases, I have found ``solver='hybr'`` and increasing the tolerance (to the same order
+as :math:`n`) help. Also be sure to update the numerical approximation values for the derivative (try something like
+``dx=10, order=15``).
+
 A warning
 -------------------
 
-Before ending this section, I want to emphasize that reducing the ``tolerance`` is not really advised. While it may
-allow the optimization routine to succeed, it only allows the 'error' of the optimization to be greater. Therefore,
-the optimization will stop 'further' from the best values. Do **not** use this approach to getting ``delicatessen`` to
-succeed in the optimization unless you are absolutely sure that the new ``tolerance`` is within acceptable limits for
-your problem. The default is `1e-9`.
+Before ending this section, I want to emphasize that reducing the ``tolerance`` is not really advised (unless you are
+in a specific scenario of a non-smooth function!). While it may allow the optimization routine to succeed, it only
+allows the 'error' of the optimization to be greater. Therefore, the optimization will stop 'further' from the best
+values. Do **not** use this approach to getting ``delicatessen`` to succeed in the optimization unless you are
+absolutely sure that the new ``tolerance`` is within acceptable limits for your problem. The default is `1e-9`.
diff --git a/docs/Reference/Estimating Equations.rst b/docs/Reference/Estimating Equations.rst
index 20a4b6f..41be6ec 100644
--- a/docs/Reference/Estimating Equations.rst
+++ b/docs/Reference/Estimating Equations.rst
@@ -28,6 +28,21 @@ Regression
 
   ee_linear_regression
   ee_logistic_regression
+  ee_robust_linear_regression
+
+
+Dose Response
+-------------------
+
+.. currentmodule:: delicatessen.estimating_equations
+
+.. autosummary::
+  :toctree: generated/
+
+  ee_4p_logistic
+  ee_3p_logistic
+  ee_2p_logistic
+  ee_effective_dose_alpha
 
 
 Causal Inference
diff --git a/docs/Reference/index.rst b/docs/Reference/index.rst
index 75d31a2..addfc73 100644
--- a/docs/Reference/index.rst
+++ b/docs/Reference/index.rst
@@ -1,7 +1,7 @@
 Reference
 ==========================
 
-This page links to the documentation for each functionality within ``delicatessen``. *M-Estimator* contains
+Documentation for all available functions and arguments for those functions are provided here. *M-Estimator* contains
 documentation for the general M-Estimator procedure. *Estimating Equations* details the built-in estimating equations
 that come with ``delicatessen``.
 
diff --git a/docs/index.rst b/docs/index.rst
index f0057bf..d47a4a0 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -3,29 +3,22 @@
 Delicatessen
 =====================================
 
-``delicatessen`` is a one-stop shop for all your sandwich (variance) needs. This Python 3.6+ package supports a
-generalization of M-Estimation theory. ``delicatessen`` draws inspiration from the R library ``geex`` for generalized
-M-estimation (`details here <https://bsaul.github.io/geex/>`_) from Saul & Hudgens (2020).
-
-Here, we provide a brief overview of M-estimation theory. For a more detailed and formal introduction to M-estimation,
-I highly recommend chapter 7 of Boos & Stefanski (2013). M-estimation is a generalization of robust inference (here
-robust refers to allowing for misspecification of secondary assumptions does not invalidate inference) for
-likelihood-based methods to a general context. *M-estimators* are solutions to estimating equations, where a generic
-estimating equation takes the form of
-
-.. math::
-
-     \sum_i^n \psi(Z_i, \theta) = 0
-
-where `Z` is independent observations (these need not be identically distributed). The parameter(s) of interest are
-the *b*-dimensional vector :math:`\theta`, and :math:`\psi` is a *b*-by-1 function that is known. A large number of
-consistent and asymptotically normal statistics can be put into the M-Estimation framework. Some examples include:
-mean, regression, delta method, and many others.
-
-To apply the M-Estimator, we solve for :math:`\theta` given the data and stacked estimating equations. This is similar
-to other approaches, but M-estimation requires the equations take the form provided above. The key advantage of
-M-Estimators is the straightforward estimation of the variance following from this framework, under suitable regularity
-conditions (and whatever else is needed). Specifically, M-Estimation provides the following sandwich variance estimator:
+``delicatessen`` is a one-stop shop for all your sandwich (variance) needs. This Python 3.6+ library supports
+M-Estimation, which is a general statistical framework for estimating unknown parameters. If you are an R user, I
+highly recommend the R library ``geex`` (`Saul & Hudgens (2020) <https://bsaul.github.io/geex/>`_).
+``delicatessen`` supports a variety of pre-built estimating equations as well as custom, user-specified estimating
+equations.
+
+Here, we provide a brief overview of M-Estimation. For a more detailed and precise introduction, please refer to
+Stefanski & Boos (2002) or Boos & Stefanski (2013). M-Estimation was developed to study the large sample properties of
+robust statistics. However, many common large-sample statistics can be expressed with estimating equations, so
+M-Estimation provides a unified structure and a streamlined approach to estimation. Let the parameter of interest be
+the vector :math:`\theta = (\theta_1, \theta_2, ..., \theta_v)` and data is observed for :math:`n` independent units
+:math:`Z_1, Z_2, …, Z_n`. Then :math:`\theta` can often be expressed as the solution to the vector equation
+:math:`\sum_{i=1}^{n} \psi(Z_i,\theta) = 0` where :math:`\psi(\dot)` is a known :math:`v \times 1`-function that does
+not depend on observation :math:`i` or :math:`n`. To compute point estimates, the vector equation is solved using the
+:math:`n` units. M-Estimators further provides a convenient and automatic method of calculating large-sample variance
+estimators via the . The sandwich variance estimator is:
 
 .. math::
 
@@ -43,14 +36,30 @@ where the prime indicates the first derivative, and
 
     B_n(Y,\hat{\theta}) = n^{-1} \sum_i^n \psi(Y_i, \hat{\theta}) \psi(Y_i, \hat{\theta})^T
 
-Therefore, we have a relatively simple way to calculate the variance for a large class of different statistics.
-Additionally, this variance estimator is robust.
+While M-Estimation is a general approach, widespread application is hindered by the corresponding derivative and matrix
+calculations. For complex estimating equations, these calculations can be especially tedious. To circumvent these
+barriers, ``delicatessen`` automates the M-Estimator.
+
+The following description is a high-level description of the process. The user provides their estimating equation(s) to
+the ``MEstimator`` class object. Next, the ``MEstimator`` object solves for :math:`\theta` using a root-finding
+algorithm. After successful completion of the root-finding step, the bread is computed by numerically approximating the
+partial derivatives and the filling is calculated. Finally, the sandwich variance is computed.
+
+
+Installation:
+-------------
+
+To install ``delicatessen``, use the following command in terminal or command prompt
+
+``python -m pip install delicatessen``
+
+Only two dependencies are necessary for ``delicatessen`` (both of which you will likely have already installed): NumPy
+and SciPy.
+
+While pandas is not necessary to use ``delicatessen``, several examples are demonstrated with pandas for ease of data
+management. To replicate the tests in ``tests/`` you will need to install ``statsmodels`` and ``pytest`` (but this is
+not necessary for general use of the package either).
 
-While M-Estimation is a powerful tool, the derivatives and matrix algebra can quickly become unwieldy. This is where
-``delicatessen`` comes in. ``delicatessen`` takes stacked estimating equations and data and works through all the necessary
-calculations. Therefore, M-Estimation can be more widely adopted without needing to solve every derivative for your
-particular problem. We can let the computer do all that hard work of finding the roots and numerically approximating the
-derivatives for us.
 
 Contents:
 -------------------------------------
@@ -65,19 +74,6 @@ Contents:
   Reference/index
   Create a GitHub Issue <https://github.com/pzivich/Deli/issues>
 
-Installation:
--------------
-
-To install ``delicatessen``, use the following command in terminal or command prompt
-
-``python -m pip install delicatessen``
-
-Only two dependencies are necessary for ``delicatessen`` (both of which you will likely have already installed): NumPy
-and SciPy.
-
-To replicate the tests in ``tests/`` you will need to install ``statsmodels`` and ``pytest`` (but this is not necessary
-for general use of the package).
-
 
 Code and Issue Tracker
 -----------------------------
diff --git a/setup.py b/setup.py
index 87bce15..7a8b11d 100644
--- a/setup.py
+++ b/setup.py
@@ -1,8 +1,9 @@
-"""Delicatessen is a one-stop shop for all your sandwich (variance) needs. Specifically, `delicatessen` implements a general API
-for M-estimation. Estimating equations are both pre-built and can be custom-made. For help on creating custom
-estimating equations, be sure to check out the ReadTheDocs documentation.
+"""Delicatessen is a one-stop shop for all your sandwich (variance) needs. `delicatessen` implements a general API for
+the generalized calculus of M-estimation. Both pre-built and custom, user-specified estimating equations are compatible
+with the M-Estimator. For an introduction to M-Estimation, using `delicatessen`, and building your own estimating
+equations, see the documentation at https://deli.readthedocs.io/en/latest/ .
 
-To install the delicatessen library, use the following command
+To install the `delicatessen`, use the following command
 
 $ python -m pip install delicatessen
 


Test patch:
diff --git a/tests/test_MEstimation.py b/tests/test_MEstimation.py
index d79e611..fb69d30 100644
--- a/tests/test_MEstimation.py
+++ b/tests/test_MEstimation.py
@@ -279,6 +279,11 @@ def psi_regression(theta):
                             np.asarray(glm.cov_params()),
                             atol=1e-6)
 
+        # Checking confidence interval estimates
+        npt.assert_allclose(mestimator.confidence_intervals(),
+                            np.asarray(glm.conf_int()),
+                            atol=1e-6)
+
     def test_logistic(self):
         """Tests linear regression by-hand with a single estimating equation.
         """
@@ -311,6 +316,11 @@ def psi_regression(theta):
                             np.asarray(glm.cov_params()),
                             atol=1e-6)
 
+        # Checking confidence interval estimates
+        npt.assert_allclose(mestimator.confidence_intervals(),
+                            np.asarray(glm.conf_int()),
+                            atol=1e-6)
+
     def test_custom_solver(self):
         """Test the use of a user-specified root-finding algorithm.
         """
diff --git a/tests/test_MEstimation_examples.py b/tests/test_MEstimation_examples.py
new file mode 100644
index 0000000..7a3d687
--- /dev/null
+++ b/tests/test_MEstimation_examples.py
@@ -0,0 +1,328 @@
+import pytest
+import numpy as np
+import pandas as pd
+from sklearn.linear_model import LinearRegression, HuberRegressor
+from delicatessen import MEstimator
+from delicatessen.estimating_equations import ee_mean, ee_mean_variance, ee_percentile, ee_positive_mean_deviation, \
+    ee_linear_regression, ee_robust_linear_regression
+
+epsilon = 1.0E-6
+
+
+class TestMEstimationExamples:
+    def test_sample_mean(self):
+        # Define dataset
+        y = np.array([-2, 1, 3, 4, 1, 4, -5, 3, 6])
+
+        def psi(theta):
+            mean = y - theta[0]
+            vari = (y - theta[0]) ** 2 - theta[1]
+            return mean, vari
+
+        mestimate = MEstimator(psi, init=[0, 1])
+        mestimate.estimate()
+
+        assert mestimate.theta[0] - np.mean(y) < epsilon
+        assert mestimate.theta[1] - np.var(y) < epsilon
+
+        # Test relative to built-in equation
+        mest = MEstimator(lambda theta: ee_mean(theta=theta, y=y), init=[0, ])
+        mest.estimate()
+
+        me = MEstimator(lambda theta: ee_mean_variance(theta, y), init=[0, 1])
+        me.estimate()
+
+        assert mestimate.theta[0] - mest.theta[0] < epsilon
+        assert mestimate.theta[0] - me.theta[0] < epsilon
+        assert mestimate.theta[1] - me.theta[1] < epsilon
+
+    # Ratio 7.2.3 --> typo in equation?
+    # Y_i - \theta_1 X_i
+    def test_ratio(self):
+        n = 10000
+        y_loc = 10
+        x_loc = 5
+        data = pd.DataFrame()
+        data['Y'] = np.random.normal(loc=y_loc, scale=2, size=n)
+        data['X'] = np.random.normal(loc=x_loc, scale=1, size=n)
+
+        def psi(theta):
+            return data['Y'] - data['X'] * theta
+
+        mestimate = MEstimator(psi, init=[1, ])
+        mestimate.estimate()
+
+        assert mestimate.theta[0] - data['Y'].mean() / data['X'].mean() < epsilon
+
+        def psi(theta):
+            mean_y = data['Y'] - theta[0]
+            mean_x = data['X'] - theta[1]
+            ratio = np.ones(data.shape[0]) * (theta[0] - theta[1] * theta[2])
+            return mean_y, mean_x, ratio
+
+        mestimate = MEstimator(psi, init=[0, 0, 1])
+        mestimate.estimate()
+
+        assert mestimate.theta[0] - data['Y'].mean() < epsilon
+        assert mestimate.theta[1] - data['X'].mean() < epsilon
+        assert mestimate.theta[2] - data['Y'].mean() / data['X'].mean() < epsilon
+
+    def test_delta_method(self):
+        n = 10000
+        y_loc = 16
+        data = pd.DataFrame(np.random.normal(loc=y_loc, scale=2, size=n),
+                            columns=['Y'])
+
+        def psi_delta(theta):
+            # Get the mean and variance (nearly identical to sample mean)
+            mean = data['Y'] - theta[0]
+            vari = (data['Y'] - theta[0]) ** 2 - theta[1]
+            # Get the standard error and log variance
+            sqrt_var = np.ones(data.shape[0]) * np.sqrt(theta[1]) - theta[2]
+            log_var = np.ones(data.shape[0]) * np.log(theta[1]) - theta[3]
+            # Return values
+            return mean, vari, sqrt_var, log_var
+
+        mestimate = MEstimator(psi_delta, init=[0, 0, 1, 1])
+        mestimate.estimate()
+
+        assert mestimate.theta[0] - data['Y'].mean() < epsilon
+        assert mestimate.theta[1] - data['Y'].var() < epsilon
+        assert mestimate.theta[2] - data['Y'].std() < epsilon
+        assert mestimate.theta[3] - np.log(data['Y'].var()) < epsilon
+
+    # Typo: Psi equations should have subscripts to index by i?
+    def test_instrumental_variable(self):
+        n = 10000
+        data = pd.DataFrame()
+        data['X'] = np.random.normal(size=n)
+        data['Y'] = 0.5 + 2 * data['X'] + np.random.normal(loc=0, size=n)
+        data['W'] = data['X'] + np.random.normal(loc=0, size=n)
+        data['T'] = -0.75 - 1 * data['X'] + np.random.normal(loc=0, size=n)
+
+        def psi(theta):
+            tinf = theta[0] - data['T']
+            tvar = (data['Y'] - theta[1] * data['W']) * (theta[0] - data['T'])
+            return tinf, tvar
+
+        def psi_instrumental(theta):
+            tinf = theta[0] - data['T']
+            winf = theta[1] - data['W']
+            wvar = (data['Y'] - theta[2] * data['W']) * (theta[1] - data['W'])
+            tvar = (data['Y'] - theta[3] * data['W']) * (theta[0] - data['T'])
+            return tinf, winf, wvar, tvar
+
+        mestimate1 = MEstimator(psi, init=[0.1, 0.1, ])
+        mestimate2 = MEstimator(psi_instrumental, init=[0.1, 0.1, 0.1, 0.1])
+        mestimate1.estimate()
+        mestimate2.estimate()
+
+        # Psi1 of docs equivalent, as is Psi2 and Psi4
+        assert mestimate1.theta[0] == mestimate2.theta[0]
+        assert mestimate1.theta[1] == mestimate2.theta[3]
+
+        # Test actual values
+
+    # Estimating equation for 7.4.1 --> provided code in example takes only lower bound, not upper bound
+    def test_robust_location(self):
+        n = 1000
+        k = 3
+        y = np.random.normal(size=n)
+
+        # Precalculate the robust mean
+        yh = np.where(y > k, k, y)
+        yhl = np.where(yh < -k, -k, yh)
+
+        def psi_robust_mean(theta):
+            k = 3
+            yr = np.where(y > k, k, y)  # Applying upper bound
+            yr = np.where(yr < -k, -k, yr)  # Applying lower bound
+            return yr - theta
+
+        estr = MEstimator(psi_robust_mean, init=[0., ])
+        estr.estimate()
+
+        print(estr.theta - yhl.mean())
+        assert estr.theta[0] - yhl.mean() < epsilon
+
+    def test_quantile_estimation(self):
+        n = 20000
+        y = np.random.normal(size=n)
+
+        # Use numpy to find quantiles
+        qs = np.quantile(y, q=[0.25, 0.50, 0.75])
+
+        def psi_quantile(theta):
+            return (
+                0.25 - 1 * (y <= theta[0]),
+                0.50 - 1 * (y <= theta[1]),
+                0.75 - 1 * (y <= theta[2])
+            )
+
+        estr = MEstimator(psi_quantile, init=[0., 0., 0., ])
+        estr.estimate(solver='hybr',
+                      tolerance=1e-4,
+                      dx=1,
+                      order=9)
+
+        # Verify using built-in equations
+        mest25 = MEstimator(lambda theta: ee_percentile(theta=theta, y=y, q=0.25), init=[-0.1, ])
+        mest50 = MEstimator(lambda theta: ee_percentile(theta=theta, y=y, q=0.50), init=[0., ])
+        mest75 = MEstimator(lambda theta: ee_percentile(theta=theta, y=y, q=0.75), init=[0.1, ])
+        mest25.estimate(solver='hybr',
+                        tolerance=1e-3,
+                        dx=1,
+                        order=15)
+        mest50.estimate(solver='hybr',
+                        tolerance=1e-3,
+                        dx=1,
+                        order=15)
+        mest75.estimate(solver='hybr',
+                        tolerance=1e-3,
+                        dx=1,
+                        order=15)
+
+        assert estr.theta[0] - mest25.theta[0] < 1e-2
+        assert estr.theta[1] - mest50.theta[0] < 1e-2
+        assert estr.theta[2] - mest75.theta[0] < 1e-2
+
+        assert estr.theta[0] - qs[0] < 1e-2
+        assert estr.theta[1] - qs[1] < 1e-2
+        assert estr.theta[2] - qs[2] < 1e-2
+
+    def test_positive_mean_deviation(self):
+        n = 10000
+        y = np.random.normal(size=n)
+
+        md = (abs(y - y.mean())).sum() / n
+
+        def psi_deviation(theta):
+            pmd = 2 * (y - theta[1]) * (y > theta[1]) - theta[0]
+            med = 0.5 - (y <= theta[1])
+            return pmd, med
+
+        estr = MEstimator(psi_deviation, init=[0., 0., ])
+        estr.estimate(solver='hybr',
+                      tolerance=1e-3,
+                      dx=1,
+                      order=9)
+
+        assert estr.theta[0] - md < 0.1
+        assert estr.theta[1] - np.median(y) < 0.1
+
+        # Test built-in equations
+        mest = MEstimator(lambda theta: ee_positive_mean_deviation(theta, y), init=[0., 0., ])
+        mest.estimate()
+
+        assert estr.theta[0] - mest.theta[0] < 0.1
+        assert estr.theta[1] - mest.theta[1] < 0.1
+
+    def test_linear_regression(self):
+        n = 1000
+        data = pd.DataFrame()
+        data['X'] = np.random.normal(size=n)
+        data['Z'] = np.random.normal(size=n)
+        data['Y'] = 0.5 + 5 * data['X'] - 2 * data['Z'] + np.random.normal(loc=0, size=n)
+        data['C'] = 1
+
+        x = np.asarray(data[['C', 'X', 'Z']])
+        y = np.asarray(data['Y'])[:, None]
+
+        model = LinearRegression().fit(x, y)
+        coef = model.coef_[0]
+
+        def psi_regression(theta):
+            beta = np.asarray(theta)[:, None]
+            return ((y - np.dot(x, beta)) * x).T
+
+        estr = MEstimator(psi_regression, init=[0., 0., 0., ])
+        estr.estimate()
+
+        assert estr.theta[0] - model.intercept_[0] < epsilon
+        assert estr.theta[1] - coef[1] < epsilon
+        assert estr.theta[2] - coef[2] < epsilon
+
+        # Test built-in equations
+        mest = MEstimator(lambda theta: ee_linear_regression(theta=theta, X=data[['C', 'X', 'Z']], y=data['Y']),
+                          init=[0., 0., 0., ])
+        mest.estimate()
+
+        assert estr.theta[0] - mest.theta[0] < epsilon
+        assert estr.theta[1] - mest.theta[1] < epsilon
+        assert estr.theta[2] - mest.theta[2] < epsilon
+
+    # Using Huber regression currently - seems to be issue with extraneous results such as
+    # extremely large estimated beta values
+    # Seems to lead to results that are equal to LR, not RobR
+    def test_robust_regression(self):
+        n = 1000
+        k = 2
+        data = pd.DataFrame()
+        data['X'] = np.random.normal(size=n)
+        data['Z'] = np.random.normal(size=n)
+        data['Y'] = 0.5 + 5 * data['X'] - 2 * data['Z'] + np.random.normal(loc=0, size=n)
+        data['C'] = 1
+
+        x = np.asarray(data[['C', 'X', 'Z']])
+        y = np.asarray(data['Y'])[:, None]
+
+        def psi_robust_regression(theta):
+            beta = np.asarray(theta)[:, None]
+            preds = np.asarray(y - np.dot(x, beta))
+            preds = np.clip(preds, a_min=-k, a_max=k)
+
+            return (preds * x).T
+
+        estr = MEstimator(psi_robust_regression, init=[0., 0., 0., ])
+        estr.estimate(solver='hybr')
+
+        model = HuberRegressor(epsilon=k).fit(x, y)
+        coef = model.coef_
+
+        assert estr.theta[0] - coef[0] < 0.3
+        assert estr.theta[1] - coef[1] < 1e-2
+        assert estr.theta[2] - coef[2] < 1e-2
+
+        # Test built-in equations
+        mest = MEstimator(lambda theta: ee_robust_linear_regression(theta, X=data[['C', 'X', 'Z']], y=data['Y'], k=2),
+                          init=[0., 0., 0., ])
+        mest.estimate(solver='hybr')
+
+        assert estr.theta[0] - mest.theta[0] < epsilon
+        assert estr.theta[1] - mest.theta[1] < epsilon
+        assert estr.theta[2] - mest.theta[2] < epsilon
+
+    # Specifically for testing whether one-lined lambda functions for built-in estimating equations throw errors.
+    def test_lambda_oneliners(self):
+        n = 10000
+        data = pd.DataFrame()
+        data['X'] = np.random.normal(size=n)
+        data['Z'] = np.random.normal(size=n)
+        data['Y'] = 0.5 + 5 * data['X'] - 2 * data['Z'] + np.random.normal(loc=0, size=n)
+        data['C'] = 1
+
+        x = np.asarray(data[['C', 'X', 'Z']])
+        y = np.asarray(data['Y'])[:, None]
+
+        # Use sklearn LinearRegression to confirm accuracy
+        model = LinearRegression().fit(x, y)
+        coef = model.coef_[0]
+
+        # Initializing MEstimator should work, but actually estimating should throw error.
+        # Guarantees that lambda function given must have variable "theta" not anything else.
+        mest = MEstimator(lambda t: ee_linear_regression(theta=t, X=data[['C', 'X', 'Z']], y=data['Y']),
+                          init=[0., 0., 0., ])
+        with pytest.raises(TypeError) as e_info:
+            mest.estimate()
+
+        # Initialize with valid lambda function
+        mest = MEstimator(lambda theta: ee_linear_regression(theta=theta, X=data[['C', 'X', 'Z']], y=data['Y']),
+                          init=[0., 0., 0., ])
+        try:
+            mest.estimate()
+        except:
+            pytest.fail("Linear regression should not throw error")
+
+        assert mest.theta[0] - model.intercept_[0] < epsilon
+        assert mest.theta[1] - coef[1] < epsilon
+        assert mest.theta[2] - coef[2] < epsilon
diff --git a/tests/test_data.py b/tests/test_data.py
new file mode 100644
index 0000000..f54681b
--- /dev/null
+++ b/tests/test_data.py
@@ -0,0 +1,16 @@
+import pytest
+
+from delicatessen.data import load_shaq_free_throws, load_inderjit
+
+
+class TestDataSets:
+
+    def test_shaq_ft_data(self):
+        d = load_shaq_free_throws()  # Loading the data
+        assert d.shape[0] == 23      # Checking number of rows
+        assert d.shape[1] == 3       # Checking number of columns
+
+    def test_inderjit(self):
+        d = load_inderjit()          # Loading the data
+        assert d.shape[0] == 24      # Checking number of rows
+        assert d.shape[1] == 2       # Checking number of columns
diff --git a/tests/test_estimating_equations.py b/tests/test_estimating_equations.py
index a134c07..4bce28b 100644
--- a/tests/test_estimating_equations.py
+++ b/tests/test_estimating_equations.py
@@ -5,15 +5,18 @@
 import statsmodels.api as sm
 import statsmodels.formula.api as smf
 from scipy.stats import logistic
+from sklearn.linear_model import Ridge
 
 from delicatessen import MEstimator
 from delicatessen.estimating_equations import (ee_mean, ee_mean_variance, ee_mean_robust,
-                                               ee_linear_regression, ee_logistic_regression,
+                                               ee_linear_regression, ee_logistic_regression, ee_ridge_linear_regression,
+                                               ee_2p_logistic, ee_3p_logistic, ee_4p_logistic, ee_effective_dose_delta,
                                                ee_gformula, ee_ipw, ee_aipw)
+from delicatessen.data import load_inderjit
 from delicatessen.utilities import inverse_logit
 
 
-class TestEstimatingEquations:
+class TestEstimatingEquationsBase:
 
     def test_mean(self):
         """Tests mean with the built-in estimating equation.
@@ -86,6 +89,9 @@ def psi2(theta):
                             mpee.asymptotic_variance,
                             atol=1e-6)
 
+
+class TestEstimatingEquationsRegression:
+
     def test_ols(self):
         """Tests linear regression with the built-in estimating equation.
         """
@@ -96,15 +102,6 @@ def test_ols(self):
         data['Y'] = 0.5 + 2*data['X'] - 1*data['Z'] + np.random.normal(loc=0, size=n)
         data['C'] = 1
 
-        def psi_regression(theta):
-            x = np.asarray(data[['C', 'X', 'Z']])
-            y = np.asarray(data['Y'])[:, None]
-            beta = np.asarray(theta)[:, None]
-            return ((y - np.dot(x, beta)) * x).T
-
-        mcee = MEstimator(psi_regression, init=[0.1, 0.1, 0.1])
-        mcee.estimate()
-
         def psi_builtin_regression(theta):
             return ee_linear_regression(theta,
                                         X=data[['C', 'X', 'Z']],
@@ -113,14 +110,22 @@ def psi_builtin_regression(theta):
         mpee = MEstimator(psi_builtin_regression, init=[0.1, 0.1, 0.1])
         mpee.estimate()
 
+        # Statsmodels function equivalent
+        glm = smf.glm("Y ~ X + Z", data).fit(cov_type="HC1")
+
         # Checking mean estimate
-        npt.assert_allclose(mcee.theta,
-                            mpee.theta,
+        npt.assert_allclose(mpee.theta,
+                            np.asarray(glm.params),
                             atol=1e-6)
 
         # Checking variance estimates
-        npt.assert_allclose(mcee.variance,
-                            mpee.variance,
+        npt.assert_allclose(mpee.variance,
+                            np.asarray(glm.cov_params()),
+                            atol=1e-6)
+
+        # Checking confidence interval estimates
+        npt.assert_allclose(mpee.confidence_intervals(),
+                            np.asarray(glm.conf_int()),
                             atol=1e-6)
 
     def test_wls(self):
@@ -151,6 +156,83 @@ def psi_regression(theta):
                             np.asarray(glm.params),
                             atol=1e-6)
 
+    def test_ridge_ols(self):
+        """Tests the ridge (L2) variation of the linear regression built-in estimating equation
+        """
+        n = 1000
+        data = pd.DataFrame()
+        data['x1'] = np.random.normal(size=n)
+        data['x2'] = data['x1'] + np.random.normal(scale=0.1, size=n)
+        data['c'] = 1
+        data['y'] = 5 + data['x1'] + np.random.normal(size=n)
+        Xvals = np.asarray(data[['c', 'x1', 'x2']])
+        yvals = np.asarray(data['y'])
+
+        # Penalty of 0.5
+        def psi(theta):
+            return ee_ridge_linear_regression(theta, X=Xvals, y=yvals, penalty=0.5, weights=None)
+
+        estr = MEstimator(psi, init=[5, 1, 1])
+        estr.estimate(solver='lm')
+        ridge_skl = Ridge(alpha=.5, fit_intercept=False).fit(X=Xvals, y=yvals)
+
+        # Checking mean estimate
+        npt.assert_allclose(estr.theta,
+                            np.asarray(ridge_skl.coef_),
+                            atol=1e-6)
+
+        # Penalty of 5.0
+        def psi(theta):
+            return ee_ridge_linear_regression(theta, X=Xvals, y=yvals, penalty=5.0, weights=None)
+
+        estr = MEstimator(psi, init=[5, 1, 1])
+        estr.estimate(solver='lm')
+        ridge_skl = Ridge(alpha=5.0, fit_intercept=False).fit(X=Xvals, y=yvals)
+
+        # Checking mean estimate
+        npt.assert_allclose(estr.theta,
+                            np.asarray(ridge_skl.coef_),
+                            atol=1e-6)
+
+    def test_ridge_wls(self):
+        """Tests the ridge (L2) variation of the weighted linear regression built-in estimating equation
+        """
+        n = 1000
+        data = pd.DataFrame()
+        data['x1'] = np.random.normal(size=n)
+        data['x2'] = data['x1'] + np.random.normal(scale=0.1, size=n)
+        data['c'] = 1
+        data['y'] = 5 + data['x1'] + np.random.normal(size=n)
+        Xvals = np.asarray(data[['c', 'x1', 'x2']])
+        yvals = np.asarray(data['y'])
+        weights = np.random.uniform(0.1, 2.5, size=n)
+
+        # Penalty of 0.5
+        def psi(theta):
+            return ee_ridge_linear_regression(theta, X=Xvals, y=yvals, penalty=0.5, weights=weights)
+
+        estr = MEstimator(psi, init=[5, 1, 1])
+        estr.estimate(solver='lm')
+        ridge_skl = Ridge(alpha=.5, fit_intercept=False).fit(X=Xvals, y=yvals, sample_weight=weights)
+
+        # Checking mean estimate
+        npt.assert_allclose(estr.theta,
+                            np.asarray(ridge_skl.coef_),
+                            atol=1e-6)
+
+        # Penalty of 5.0
+        def psi(theta):
+            return ee_ridge_linear_regression(theta, X=Xvals, y=yvals, penalty=5.0, weights=weights)
+
+        estr = MEstimator(psi, init=[5, 1, 1])
+        estr.estimate(solver='lm')
+        ridge_skl = Ridge(alpha=5.0, fit_intercept=False).fit(X=Xvals, y=yvals, sample_weight=weights)
+
+        # Checking mean estimate
+        npt.assert_allclose(estr.theta,
+                            np.asarray(ridge_skl.coef_),
+                            atol=1e-6)
+
     def test_logitic(self):
         n = 500
         data = pd.DataFrame()
@@ -159,15 +241,6 @@ def test_logitic(self):
         data['Y'] = np.random.binomial(n=1, p=logistic.cdf(0.5 + 2*data['X'] - 1*data['Z']), size=n)
         data['C'] = 1
 
-        def psi_regression(theta):
-            x = np.asarray(data[['C', 'X', 'Z']])
-            y = np.asarray(data['Y'])[:, None]
-            beta = np.asarray(theta)[:, None]
-            return ((y - inverse_logit(np.dot(x, beta))) * x).T
-
-        mcee = MEstimator(psi_regression, init=[0., 0., 0.])
-        mcee.estimate()
-
         def psi_builtin_regression(theta):
             return ee_logistic_regression(theta,
                                           X=data[['C', 'X', 'Z']],
@@ -176,14 +249,22 @@ def psi_builtin_regression(theta):
         mpee = MEstimator(psi_builtin_regression, init=[0., 0., 0.])
         mpee.estimate()
 
+        # Comparing to statsmodels GLM (with robust covariance)
+        glm = smf.glm("Y ~ X + Z", data, family=sm.families.Binomial()).fit(cov_type="HC1")
+
         # Checking mean estimate
-        npt.assert_allclose(mcee.theta,
-                            mpee.theta,
+        npt.assert_allclose(mpee.theta,
+                            np.asarray(glm.params),
                             atol=1e-6)
 
         # Checking variance estimates
-        npt.assert_allclose(mcee.variance,
-                            mpee.variance,
+        npt.assert_allclose(mpee.variance,
+                            np.asarray(glm.cov_params()),
+                            atol=1e-6)
+
+        # Checking confidence interval estimates
+        npt.assert_allclose(mpee.confidence_intervals(),
+                            np.asarray(glm.conf_int()),
                             atol=1e-6)
 
     def test_weighted_logistic(self):
@@ -214,6 +295,180 @@ def psi_regression(theta):
                             np.asarray(glm.params),
                             atol=1e-6)
 
+
+class TestEstimatingEquationsDoseResponse:
+
+    def test_4pl(self):
+        """Test the 4 parameter log-logistic model using Inderjit et al. (2002)
+
+        Compares against R's drc library:
+
+        library(drc)
+        library(sandwich)
+        library(lmtest)
+
+        data(ryegrass)
+        rgll4 = drm(rootl ~ conc, data=ryegrass, fct=LL.4())
+        coeftest(rgll4, vcov=sandwich)
+        """
+        d = load_inderjit()
+        dose_data = d[:, 1]
+        resp_data = d[:, 0]
+
+        def psi(theta):
+            return ee_4p_logistic(theta=theta, X=dose_data, y=resp_data)
+
+        # Optimization procedure
+        mestimator = MEstimator(psi, init=[0, 2, 1, 10])
+        mestimator.estimate(solver='lm')
+
+        # R optimization from Ritz et al.
+        comparison_theta = np.asarray([0.48141, 3.05795, 2.98222, 7.79296])
+        comparison_var = np.asarray([0.12779, 0.26741, 0.47438, 0.15311])
+
+        # Checking mean estimate
+        npt.assert_allclose(mestimator.theta,
+                            comparison_theta,
+                            atol=1e-5)
+
+        # Checking variance estimate
+        npt.assert_allclose(np.diag(mestimator.variance)**0.5,
+                            comparison_var,
+                            atol=1e-4)
+
+    def test_3pl(self):
+        """Test the 3 parameter log-logistic model using Inderjit et al. (2002)
+
+        Compares against R's drc library:
+
+        library(drc)
+        library(sandwich)
+        library(lmtest)
+
+        data(ryegrass)
+        rgll3 = drm(rootl ~ conc, data=ryegrass, fct=LL.3())
+        coeftest(rgll3, vcov=sandwich)
+        """
+        d = load_inderjit()
+        dose_data = d[:, 1]
+        resp_data = d[:, 0]
+
+        def psi(theta):
+            return ee_3p_logistic(theta=theta, X=dose_data, y=resp_data,
+                                  lower=0)
+
+        # Optimization procedure
+        mestimator = MEstimator(psi, init=[2, 1, 10])
+        mestimator.estimate(solver='lm')
+
+        # R optimization from Ritz et al.
+        comparison_theta = np.asarray([3.26336, 2.47033, 7.85543])
+        comparison_var = np.asarray([0.26572, 0.29238, 0.15397])
+
+        # Checking mean estimate
+        npt.assert_allclose(mestimator.theta,
+                            comparison_theta,
+                            atol=1e-5)
+
+        # Checking variance estimate
+        npt.assert_allclose(np.diag(mestimator.variance)**0.5,
+                            comparison_var,
+                            atol=1e-5)
+
+    def test_2pl(self):
+        """Test the 2 parameter log-logistic model using Inderjit et al. (2002)
+
+        Compares against R's drc library:
+
+        library(drc)
+        library(sandwich)
+        library(lmtest)
+
+        data(ryegrass)
+        rgll2 = drm(rootl ~ conc, data=ryegrass, fct=LL.2(upper=8))
+        coeftest(rgll2, vcov=sandwich)
+        """
+        d = load_inderjit()
+        dose_data = d[:, 1]
+        resp_data = d[:, 0]
+
+        def psi(theta):
+            return ee_2p_logistic(theta=theta, X=dose_data, y=resp_data,
+                                  lower=0, upper=8)
+
+        # Optimization procedure
+        mestimator = MEstimator(psi, init=[2, 1])
+        mestimator.estimate(solver='lm')
+
+        # R optimization from Ritz et al.
+        comparison_theta = np.asarray([3.19946, 2.38220])
+        comparison_var = np.asarray([0.24290, 0.27937])
+
+        # Checking mean estimate
+        npt.assert_allclose(mestimator.theta,
+                            comparison_theta,
+                            atol=1e-5)
+
+        # Checking variance estimate
+        npt.assert_allclose(np.diag(mestimator.variance)**0.5,
+                            comparison_var,
+                            atol=1e-5)
+
+    def test_3pl_ed_delta(self):
+        """Test the ED(alpha) calculation with the 3 parameter log-logistic model using Inderjit et al. (2002)
+
+        Compares against R's drc library:
+
+        library(drc)
+        library(sandwich)
+
+        data(ryegrass)
+        rgll3 = drm(rootl ~ conc, data=ryegrass, fct=LL.3())
+        ED(rgll3, c(5, 10, 50), interval='delta', vcov=sandwich)
+        """
+        d = load_inderjit()
+        dose_data = d[:, 1]
+        resp_data = d[:, 0]
+
+        def psi(theta):
+            lower_limit = 0
+            pl3 = ee_3p_logistic(theta=theta, X=dose_data, y=resp_data,
+                                 lower=lower_limit)
+            ed05 = ee_effective_dose_delta(theta[3], y=resp_data, delta=0.05,
+                                           steepness=theta[0], ed50=theta[1],
+                                           lower=lower_limit, upper=theta[2])
+            ed10 = ee_effective_dose_delta(theta[4], y=resp_data, delta=0.10,
+                                           steepness=theta[0], ed50=theta[1],
+                                           lower=lower_limit, upper=theta[2])
+            ed50 = ee_effective_dose_delta(theta[5], y=resp_data, delta=0.50,
+                                           steepness=theta[0], ed50=theta[1],
+                                           lower=lower_limit, upper=theta[2])
+            return np.vstack((pl3,
+                              ed05,
+                              ed10,
+                              ed50))
+
+        # Optimization procedure
+        mestimator = MEstimator(psi, init=[2, 1, 10, 1, 1, 2])
+        mestimator.estimate(solver='lm')
+
+        # R optimization from Ritz et al.
+        comparison_theta = np.asarray([0.99088, 1.34086, 3.26336])
+        comparison_var = np.asarray([0.12397, 0.13134, 0.26572])
+
+        # Checking mean estimate
+        npt.assert_allclose(mestimator.theta[-3:],
+                            comparison_theta,
+                            atol=1e-5)
+
+        # Checking variance estimate
+        npt.assert_allclose(np.diag(mestimator.variance)[-3:]**0.5,
+                            comparison_var,
+                            atol=1e-5)
+
+
+class TestEstimatingEquationsCausal:
+
     @pytest.fixture
     def causal_data(self):
         np.random.seed(1205811)
@@ -233,12 +488,18 @@ def causal_data(self):
         return df
 
     def test_gformula(self, causal_data):
+        d1 = causal_data.copy()
+        d1['A'] = 1
+        d0 = causal_data.copy()
+        d0['A'] = 0
+
         # M-estimation
         def psi(theta):
             return ee_gformula(theta,
-                               X=causal_data[['C', 'A', 'W']],
                                y=causal_data['Y'],
-                               treat_index=1)
+                               X=causal_data[['C', 'A', 'W']],
+                               X1=d1[['C', 'A', 'W']],
+                               X0=d0[['C', 'A', 'W']])
 
         mestimator = MEstimator(psi, init=[0., 0.5, 0.5, 0., 0., 0.])
         mestimator.estimate(solver='lm')
@@ -267,13 +528,75 @@ def psi(theta):
                             np.mean(ya0),
                             atol=1e-6)
 
+    def test_gcomp_bad_dimensions_error(self, causal_data):
+        d1 = causal_data.copy()
+        d1['A'] = 1
+        d0 = causal_data.copy()
+        d0['A'] = 0
+
+        # M-estimation
+        def psi(theta):
+            return ee_gformula(theta,
+                               y=causal_data['Y'],
+                               X=causal_data[['C', 'A', 'W']],
+                               X1=d1[['C', 'W']])
+
+        mestimator = MEstimator(psi, init=[0.5, 0., 0., 0.])
+        with pytest.raises(ValueError, match="The dimensions of X and X1"):
+            mestimator.estimate(solver='lm')
+
+        def psi(theta):
+            return ee_gformula(theta,
+                               y=causal_data['Y'],
+                               X=causal_data[['C', 'A', 'W']],
+                               X1=d1[['C', 'A', 'W']],
+                               X0=d0[['C', 'A']])
+
+        mestimator = MEstimator(psi, init=[0., 0.5, 0.5, 0., 0., 0.])
+        with pytest.raises(ValueError, match="The dimensions of X and X0"):
+            mestimator.estimate(solver='lm')
+
     def test_ipw(self, causal_data):
         # M-estimation
         def psi(theta):
             return ee_ipw(theta,
-                          X=causal_data[['C', 'A', 'W']],
                           y=causal_data['Y'],
-                          treat_index=1)
+                          A=causal_data['A'],
+                          W=causal_data[['C', 'W']])
+
+        mestimator = MEstimator(psi, init=[0., 0.5, 0.5, 0., 0.])
+        mestimator.estimate(solver='lm')
+
+        # By-hand IPW estimator with statsmodels
+        glm = sm.GLM(causal_data['A'], causal_data[['C', 'W']],
+                     family=sm.families.Binomial()).fit()
+        pi = glm.predict()
+        ya1 = causal_data['A'] * causal_data['Y'] / pi
+        ya0 = (1-causal_data['A']) * causal_data['Y'] / (1-pi)
+
+        # Checking logistic coefficients (nuisance model estimates)
+        npt.assert_allclose(mestimator.theta[3:],
+                            np.asarray(glm.params),
+                            atol=1e-6)
+        # Checking mean estimates
+        npt.assert_allclose(mestimator.theta[0],
+                            np.mean(ya1) - np.mean(ya0),
+                            atol=1e-6)
+        npt.assert_allclose(mestimator.theta[1],
+                            np.mean(ya1),
+                            atol=1e-6)
+        npt.assert_allclose(mestimator.theta[2],
+                            np.mean(ya0),
+                            atol=1e-6)
+
+    def test_ipw_truncate(self, causal_data):
+        # M-estimation
+        def psi(theta):
+            return ee_ipw(theta,
+                          y=causal_data['Y'],
+                          A=causal_data['A'],
+                          W=causal_data[['C', 'W']],
+                          truncate=(0.1, 0.5))
 
         mestimator = MEstimator(psi, init=[0., 0.5, 0.5, 0., 0.])
         mestimator.estimate(solver='lm')
@@ -282,6 +605,7 @@ def psi(theta):
         glm = sm.GLM(causal_data['A'], causal_data[['C', 'W']],
                      family=sm.families.Binomial()).fit()
         pi = glm.predict()
+        pi = np.clip(pi, 0.1, 0.5)
         ya1 = causal_data['A'] * causal_data['Y'] / pi
         ya0 = (1-causal_data['A']) * causal_data['Y'] / (1-pi)
 
@@ -300,13 +624,34 @@ def psi(theta):
                             np.mean(ya0),
                             atol=1e-6)
 
+    def test_ipw_truncate_error(self, causal_data):
+        # M-estimation
+        def psi(theta):
+            return ee_ipw(theta,
+                          y=causal_data['Y'],
+                          A=causal_data['A'],
+                          W=causal_data[['C', 'W']],
+                          truncate=(0.99, 0.01))
+
+        mestimator = MEstimator(psi, init=[0., 0.5, 0.5, 0., 0.])
+        with pytest.raises(ValueError, match="truncate values"):
+            mestimator.estimate()
+
     def test_aipw(self, causal_data):
+        d1 = causal_data.copy()
+        d1['A'] = 1
+        d0 = causal_data.copy()
+        d0['A'] = 0
+
         # M-estimation
         def psi_builtin_regression(theta):
             return ee_aipw(theta,
-                           X=causal_data[['C', 'A', 'W']],
                            y=causal_data['Y'],
-                           treat_index=1)
+                           A=causal_data['A'],
+                           W=causal_data[['C', 'W']],
+                           X=causal_data[['C', 'A', 'W']],
+                           X1=d1[['C', 'A', 'W']],
+                           X0=d0[['C', 'A', 'W']])
 
         mestimator = MEstimator(psi_builtin_regression, init=[0., 0.5, 0.5,   # Parameters of interest
                                                               0., 0., 0.,     # Outcome nuisance model
@@ -334,12 +679,12 @@ def psi_builtin_regression(theta):
         var_r0 = np.nanvar(ya0_star - np.mean(ya0_star), ddof=1) / causal_data.shape[0]
 
         # Checking logistic coefficients (nuisance model estimates)
-        npt.assert_allclose(mestimator.theta[3:6],
-                            np.asarray(y_m.params),
-                            atol=1e-6)
-        npt.assert_allclose(mestimator.theta[6:],
+        npt.assert_allclose(mestimator.theta[3:5],
                             np.asarray(pi_m.params),
                             atol=1e-6)
+        npt.assert_allclose(mestimator.theta[5:],
+                            np.asarray(y_m.params),
+                            atol=1e-6)
 
         # Checking mean estimates
         npt.assert_allclose(mestimator.theta[0],